Thank you very much! I’ll do it hoping that it will be helpful.
Do not set “best” cvx precision. Use default precision and everything solves nicely (I’m talking about Mosek).
As discussed in some other thread I can’t find any more, highest cvx precision sets Mosek termination tolerances to 0.0, which strictly speaking should be illegal and can yield all kinds of crazy behavior.
Ok, thank you, I will remember it from now on.
The problem is that now the solution is the same given by the other solvers, and for some odd reason it violates one of the two constraints of the problem.
Hello to everyone,
I’m now pretty sure I found out where the problem is.
Basically, Inside of my cvx workspace I declare two variables:
cvx_begin %quiet
% Design variable declaration
variable K1(N*nu,nx);
variable Ktri(N*nu,nx,N+1);
which dimensions are not important.
I use K1 and Ktri to create a variable K and some other cvx expressions.
Now, as @Mark_L_Stone told me, I’m not guaranteed to get the right expressions after cvx, is closed, but I should compute everything back starting from K1 and Ktri.
The odd fact is that if I do it, the values I get are wrong ( just the most critical example, the value of the cost function computed with the output K1 and Ktri is lower than cvx_optval, and the initial issue (not satisfying the constraint)was just a consequence).
I notice that because adding a very small term to the cost function (of the order of 10^-100) I got the same cvx_optval but the output K1 and Ktri were such that the constraint was satisfied.
Therefore my question now becomes:
How is it possible that the output design variables are not the right ones and how can I be sure to get the right ones?
Thank you in advance.
CVX variables should have the correct values.If this does not appear to be the case, please show us the detailed information leading you to that conclusion.
Hello, here is the CVX problem and the code I run right after to compare the optimal value computed by CVX and the one obtained from the CVX variable given as output :
cvx_begin %quiet
variable K1(N*nu,nx);
variable Ktri(N*nu,nx,N+1);
KK=[K1];
for ii=1:N+1
KK=[KK blkdiag(zeros(ii-1)*nu,eye((N-(ii-1))*nu))*Ktri(:,:,ii)];
end
K=KK;
%% Cost function
% Cost function elements for compostion
U=mu_0_aug*mu_0_aug';
H1=AAA*U*AAA';
V=DDD*DDD';
Oo=EEE*Sig_0_aug*EEE';
Z=V+Oo;
n1=trace(QQQ*Z);
n2=trace(K'*BBB'*QQQ*Z);
n3=trace(QQQ*BBB*K*Z);
soft1 = trace(SSSc*H1);
soft2 = trace(K'*BBB'*SSSc*H1);
soft3 = trace(SSSc*BBB*K*H1);
n4=2*trace(QQQ*((I+BBB*K)*(EEE*Sig_0_aug*AAA_errev'+DDD*DDD_errev'))); %cross covariance term
tmp6=0;
F6=BBB*K*DDD;
for k=1:size(F6,2)
tmp6=tmp6+quad_form(F6(:,k),QQQ);
end
n6=tmp6;
tmp7=0;
F7=BBB*K*EEE*(Sig_0_aug)^0.5;
for k=1:size(F7,2)
tmp7=tmp7+quad_form(F7(:,k),QQQ);
end
n7=tmp7;
tmp8=0;
F8=K*AAA*mu_0_aug;
for k=1:size(F8,2)
tmp8=tmp8+quad_form(F8(:,k),RRR);
end
n8=tmp8;
tmp10=0;
F10=K*DDD;
for k=1:size(F10,2)
tmp10=tmp10+quad_form(F10(:,k),RRR);
end
n10=tmp10;
tmp11=0;
F11=K*EEE*(Sig_0_aug)^0.5;
for k=1:size(F11,2)
tmp11=tmp11+quad_form(F11(:,k),RRR);
end
n11=tmp11;
tmp_s4=0;
Fs4=BBB*K*AAA*mu_0_aug;
for k=1:size(Fs4,2)
tmp_s4=tmp_s4+quad_form(Fs4(:,k),SSSc);
end
soft4=tmp_s4;
%Cost function
J=n1+n2+n3+n4+n6+n7+n8+n10+n11+soft1+soft2+soft3+soft4;
minimize(J);
subject to
beta-alpha'*(I+BBB*K)*AAA*mu_0_aug-(norm(((Sig_op)^0.5)*(I+BBB*K)'*alpha,2)+norm((Sig_aug)^0.5*alpha,2))*failnorm <= 0;
1-norm( (Sig_op^0.5*(I+BBB*K)'*E_N'*Sig_lim^(-0.5)) ,2)>=0;
cvx_end
%Rebuilding of the cost function after the optimization is over
KK=[K1];
for ii=1:N+1
KK=[KK blkdiag(zeros(ii-1)*nu,eye((N-(ii-1))*nu))*Ktri(:,:,ii)];
end
K=KK;
%% Cost function
% Cost function elements for compostion
U=mu_0_aug*mu_0_aug';
H1=AAA*U*AAA';
V=DDD*DDD';
Oo=EEE*Sig_0_aug*EEE';
Z=V+Oo;
n1=trace(QQQ*Z);
n2=trace(K'*BBB'*QQQ*Z);
n3=trace(QQQ*BBB*K*Z);
soft1 = trace(SSSc*H1);
soft2 = trace(K'*BBB'*SSSc*H1);
soft3 = trace(SSSc*BBB*K*H1);
n4=2*trace(QQQ*((I+BBB*K)*(EEE*Sig_0_aug*AAA_errev'+DDD*DDD_errev'))); %cross covariance term
tmp6=0;
F6=BBB*K*DDD;
for k=1:size(F6,2)
tmp6=tmp6+quad_form(F6(:,k),QQQ);
end
n6=tmp6;
tmp7=0;
F7=BBB*K*EEE*(Sig_0_aug)^0.5;
for k=1:size(F7,2)
tmp7=tmp7+quad_form(F7(:,k),QQQ);
end
n7=tmp7;
tmp8=0;
F8=K*AAA*mu_0_aug;
for k=1:size(F8,2)
tmp8=tmp8+quad_form(F8(:,k),RRR);
end
n8=tmp8;
tmp10=0;
F10=K*DDD;
for k=1:size(F10,2)
tmp10=tmp10+quad_form(F10(:,k),RRR);
end
n10=tmp10;
tmp11=0;
F11=K*EEE*(Sig_0_aug)^0.5;
for k=1:size(F11,2)
tmp11=tmp11+quad_form(F11(:,k),RRR);
end
n11=tmp11;
tmp_s4=0;
Fs4=BBB*K*AAA*mu_0_aug;
for k=1:size(Fs4,2)
tmp_s4=tmp_s4+quad_form(Fs4(:,k),SSSc);
end
soft4=tmp_s4;
%Cost function
J=n1+n2+n3+n4+n6+n7+n8+n10+n11+soft1+soft2+soft3+soft4;
As you can see I rewrite the exact code that was internal to the CVX environment starting from the CVX variables K1 and Ktri.
The log output from CVX (using Mosek) is:
Calling Mosek 9.1.9: 12679 variables, 827 equality constraints
For improved efficiency, Mosek is solving the dual problem.
------------------------------------------------------------
MOSEK Version 9.1.9 (Build date: 2019-11-21 11:34:40)
Copyright (c) MOSEK ApS, Denmark. WWW: mosek.com
Platform: Windows/64-X86
MOSEK warning 710: #1 (nearly) zero elements are specified in sparse col '' (15) of matrix 'A'.
MOSEK warning 710: #1 (nearly) zero elements are specified in sparse col '' (17) of matrix 'A'.
MOSEK warning 710: #3 (nearly) zero elements are specified in sparse col '' (18) of matrix 'A'.
MOSEK warning 710: #3 (nearly) zero elements are specified in sparse col '' (20) of matrix 'A'.
MOSEK warning 710: #6 (nearly) zero elements are specified in sparse col '' (21) of matrix 'A'.
MOSEK warning 710: #1 (nearly) zero elements are specified in sparse col '' (22) of matrix 'A'.
MOSEK warning 710: #6 (nearly) zero elements are specified in sparse col '' (23) of matrix 'A'.
MOSEK warning 710: #10 (nearly) zero elements are specified in sparse col '' (24) of matrix 'A'.
MOSEK warning 710: #4 (nearly) zero elements are specified in sparse col '' (25) of matrix 'A'.
MOSEK warning 710: #11 (nearly) zero elements are specified in sparse col '' (26) of matrix 'A'.
Warning number 710 is disabled.
Problem
Name :
Objective sense : min
Type : CONIC (conic optimization problem)
Constraints : 827
Cones : 236
Scalar variables : 10264
Matrix variables : 1
Integer variables : 0
Optimizer started.
Presolve started.
Linear dependency checker started.
Linear dependency checker terminated.
Eliminator started.
Freed constraints in eliminator : 1
Eliminator terminated.
Eliminator started.
Freed constraints in eliminator : 0
Eliminator terminated.
Eliminator - tries : 2 time : 0.00
Lin. dep. - tries : 1 time : 0.00
Lin. dep. - number : 0
Presolve terminated. Time: 0.08
Problem
Name :
Objective sense : min
Type : CONIC (conic optimization problem)
Constraints : 827
Cones : 236
Scalar variables : 10264
Matrix variables : 1
Integer variables : 0
Optimizer - threads : 4
Optimizer - solved problem : the primal
Optimizer - Constraints : 807
Optimizer - Cones : 236
Optimizer - Scalar variables : 4976 conic : 4975
Optimizer - Semi-definite variables: 1 scalarized : 2415
Factor - setup time : 0.03 dense det. time : 0.00
Factor - ML order time : 0.00 GP order time : 0.00
Factor - nonzeros before factor : 2.10e+05 after factor : 2.13e+05
Factor - dense dim. : 0 flops : 2.38e+08
ITE PFEAS DFEAS GFEAS PRSTATUS POBJ DOBJ MU TIME
0 6.9e+01 4.0e+00 2.0e+00 0.00e+00 -3.043578051e+00 0.000000000e+00 1.0e+00 0.19
1 2.9e+01 1.7e+00 1.1e+00 -7.34e-01 -6.855128468e+01 -6.535390748e+01 4.2e-01 0.36
2 1.7e+01 9.7e-01 5.7e-01 -2.55e-01 -2.171738190e+02 -2.145765599e+02 2.4e-01 0.42
3 7.6e+00 4.5e-01 2.1e-01 1.66e-01 -4.987466018e+02 -4.972026950e+02 1.1e-01 0.50
4 2.2e+00 1.3e-01 3.3e-02 3.56e-01 -7.670074253e+02 -7.665184783e+02 3.2e-02 0.59
5 1.2e-01 7.1e-03 6.0e-04 7.67e-01 -9.346072680e+02 -9.345639833e+02 1.7e-03 0.67
6 1.4e-03 8.4e-05 8.1e-07 9.87e-01 -9.449396656e+02 -9.449391249e+02 2.1e-05 0.77
7 1.4e-05 8.1e-07 7.6e-10 1.00e+00 -9.450620293e+02 -9.450620241e+02 2.0e-07 0.84
8 1.3e-07 7.4e-09 6.5e-13 1.00e+00 -9.450632083e+02 -9.450632082e+02 1.8e-09 0.91
Optimizer terminated. Time: 0.98
Interior-point solution summary
Problem status : PRIMAL_AND_DUAL_FEASIBLE
Solution status : OPTIMAL
Primal. obj: -9.4506320829e+02 nrm: 1e+02 Viol. con: 5e-07 var: 0e+00 barvar: 0e+00 cones: 0e+00
Dual. obj: -9.4506320824e+02 nrm: 1e+01 Viol. con: 0e+00 var: 8e-09 barvar: 8e-09 cones: 0e+00
Optimizer summary
Optimizer - time: 0.98
Interior-point - iterations : 8 time: 0.91
Basis identification - time: 0.00
Primal - iterations : 0 time: 0.00
Dual - iterations : 0 time: 0.00
Clean primal - iterations : 0 time: 0.00
Clean dual - iterations : 0 time: 0.00
Simplex - time: 0.00
Primal simplex - iterations : 0 time: 0.00
Dual simplex - iterations : 0 time: 0.00
Mixed integer - relaxations: 0 time: 0.00
------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): +162.303
where you can see that cvx_optval =162.303, while the variable J that is computed later is equal to 117.9348.
If I slightly change one matrix inside the problem before the optimization (Sig_op, which elements are of the order of 10^(-2)), adding to it a matrix that I think is negligible (its biggest element is of the order of 10^(-98) ) and I run the same code, I get the same cvx_optval =162.303, but I get also J=162.3030 as it should be.
Did you send the data to Mosek that they suggested?
I warned you before about the Mosek warnings. Perhaps those near zero elements are causing numerical inaccuracy and sensitivity.
I didn’t verify your reconstruction of the optimal objective value is correct, but I’m guessing so. Unless there is a bug somewhere, it appears that numerical difficulties due to “bad” input data is to blame for your observed discrepancy.
The Mosek solution summary says
Interior-point solution summary
Problem status : PRIMAL_AND_DUAL_FEASIBLE
Solution status : OPTIMAL
Primal. obj: -9.4506320829e+02 nrm: 1e+02 Viol. con: 5e-07 var: 0e+00 barvar: 0e+00 cones: 0e+00
Dual. obj: -9.4506320824e+02 nrm: 1e+01 Viol. con: 0e+00 var: 8e-09 barvar: 8e-09 cones: 0e+00
that means Mosek solves the problem accurately. I am 99.99% sure the problem is not with Mosek.
Your problem or Cvx is culprit.
Btw how come 10^(-98) pops up in your argument i.e.
(its biggest element is of the order of 10^(-98) )
? You are aware
1+epsilon=1
on a computer for |epsilon|<=1.0e-16. So mentioning 10^(-98) seems very strange, In my numerical analysis course I learned using 10^(-98) was most likely a very bad idea.