Dual Infeasibility Certificate and Bounding SDPs
216 views
Skip to first unread message
Chris D'Angelo
Apr 11, 2018, 9:49:43 AM4/11/18
to YALMIP
Dr. Löfberg,
I have a modeling question that I believe you can help me with. I am using YALMIP with MOSEK to solve SDPs. Specifically, I am synthesizing output-feedback H-infinity controllers for uncertain plants. In the first portion of my code, which I will not include here since it is not causing me any problems, I generate random instances of my plant model and solve for two PSD variables, which I will call R and S, that indicate the existence of an H-infinity controller. To achieve this end, additional LMI constraints, representative of instances of this random plant, are appended as constraints in this SDP as a feasibility problem.
I am able to find a feasible R and S for my uncertain system. From here, I must construct the controller from these matrices (Gahinet/Apkirian 1994).
%%%%%%%%%%%%%%
k = rank(eye(n) - Ropt*Sopt);
R = Ropt;
S = Sopt;
Theta = sdpvar(k+m2,k+p2,'full');
G = [];
for r = 1:Scenarios
%%%% GET SOME MATRIX SIZES %%%%
n = size(A,2);
m1 = size(B1,2);
m2 = size(B2,2);
p1 = size(C1,1);
p2 = size(C2,1);
[U,Sv,Vt] = svd(eye(n) - R*S); %Equation 7.1
M = U; Nt = (Sv*Vt');
%Now, we move into solving equation 7.2
%We note that Xcl is (n+k)x(n+k).
X1 = [R eye(size(R,1));...
M' zeros(size(M',1),size(R,2))];
X2 = [eye(size(S,1)) S;...
zeros(size(Nt,1),size(S,1)) Nt];
Xcl = X2/X1;
%Now, we need to formulate 7.3 and solve another LMI to get our controller
%variables
A0 = [As(:,:,r) zeros(n,k);...
zeros(k,n) zeros(k)];
B0 = [B1s(:,:,r);zeros(k,size(B1,2))];
C0 = [C1s(:,:,r) zeros(size(C1,1),k)];
Bcal = [zeros(size(B2,1),k) B2s(:,:,r);...
eye(k) zeros(k,size(B2,2))];
Ccal = [zeros(k,size(C2,2)),eye(k);...
C2s(:,:,r), zeros(size(C2,1),k)];
D12cal = [zeros(p1,k),D12];
D21cal = [zeros(k,m1);D21];
%These are per equations in 4.2
PcalXcl = [Bcal'*Xcl zeros(k+m2,m1) D12cal']; %From equation 4.7
Qcal = [Ccal D21cal zeros(k+p2,p1)];
%This is from equation 4.4
PSI_Xcl = [A0'*Xcl+Xcl*A0 Xcl*B0 C0';...
B0'*Xcl -gamma*eye(size(Xcl*B0,2)) D11';...
C0 D11 -gamma*eye(size(D11,1))];
Fcontroller = PSI_Xcl + Qcal'*Theta'*PcalXcl + PcalXcl'*Theta*Qcal;
Fcontroller = 0.5*(Fcontroller+Fcontroller')<=0;
G = [G Fcontroller];
end
ops = sdpsettings('solver','mosek','verbose',1,'mosek.MSK_IPAR_INFEAS_REPORT_AUTO','MSK_ON','mosek.MSK_IPAR_INFEAS_REPORT_LEVEL',1,'savesolveroutput',1,'dualize',0);
sol_controller = optimize(G,0,ops);
%%%%%%%%%%%%%%
The code snippet shown here is not executable, but I wanted you to see how the sdpvar was called and how the matrix variable in the problem was formulated. Ultimately, Theta = [Ak Bk; Ck Dk], which are the state space versions of my controller matrices. Again, R and S were found such that they indicated the existence of a controller for the random plant scenarios indexed by As(:,:,:), Bs1(:,:,:), etc. Therefore, I am attempting to construct the controller matrices, contained within Theta, by solving this SDP. Here is the output that I get (for Scenarios = 10):
Problem
Name :
Objective sense : min
Type : CONIC (conic optimization problem)
Constraints : 504
Cones : 0
Scalar variables : 0
Matrix variables : 10
Integer variables : 0
Optimizer started.
Presolve started.
Linear dependency checker started.
Linear dependency checker terminated.
Eliminator - tries : 0 time : 0.00
Lin. dep. - tries : 1 time : 0.00
Lin. dep. - number : 0
Presolve terminated. Time: 0.00
Problem
Name :
Objective sense : min
Type : CONIC (conic optimization problem)
Constraints : 504
Cones : 0
Scalar variables : 0
Matrix variables : 10
Integer variables : 0
Optimizer - threads : 8
Optimizer - solved problem : the primal
Optimizer - Constraints : 504
Optimizer - Cones : 0
Optimizer - Scalar variables : 0 conic : 0
Optimizer - Semi-definite variables: 10 scalarized : 12750
Factor - setup time : 0.04 dense det. time : 0.00
Factor - ML order time : 0.01 GP order time : 0.00
Factor - nonzeros before factor : 1.27e+05 after factor : 1.27e+05
Factor - dense dim. : 0 flops : 5.48e+08
ITE PFEAS DFEAS GFEAS PRSTATUS POBJ DOBJ MU TIME
0 1.0e+01 1.2e+04 3.3e+01 0.00e+00 3.313359611e+04 0.000000000e+00 1.0e+00 0.10
1 8.5e-02 9.9e+01 2.6e-02 -1.00e+00 -8.623959936e+04 0.000000000e+00 8.5e-03 0.37
2 1.7e-02 1.9e+01 2.4e-03 -9.96e-01 -5.238302674e+05 0.000000000e+00 1.7e-03 0.55
3 1.3e-03 1.5e+00 5.0e-05 -9.99e-01 -7.357798193e+06 0.000000000e+00 1.3e-04 0.74
4 1.0e-04 1.2e-01 1.2e-06 -1.00e+00 -9.000527122e+07 0.000000000e+00 1.0e-05 0.93
5 2.6e-05 3.1e-02 1.5e-07 -1.00e+00 -3.501866241e+08 0.000000000e+00 2.6e-06 1.11
6 4.5e-06 5.3e-03 1.1e-08 -1.00e+00 -2.056547747e+09 0.000000000e+00 4.5e-07 1.30
7 6.3e-08 7.3e-05 1.7e-11 -1.00e+00 -1.482900559e+11 0.000000000e+00 6.3e-09 1.49
8 3.4e-12 7.8e-16 7.7e-10 -1.00e+00 -9.341766352e+02 0.000000000e+00 4.9e-17 1.67
Optimizer terminated. Time: 1.71
MOSEK DUAL INFEASIBILITY REPORT.
Problem status: The problem is dual infeasible
The following constraints are involved in the infeasibility.
Index Name Activity Objective Lower bound Upper bound
The following variables are involved in the infeasibility.
Index Name Activity Objective Lower bound Upper bound
Interior-point solution summary
Problem status : DUAL_INFEASIBLE
Solution status : DUAL_INFEASIBLE_CER
Primal. obj: -9.3417663519e+02 nrm: 2e+02 Viol. con: 9e-09 barvar: 0e+00
Optimizer summary
Optimizer - time: 1.71
Interior-point - iterations : 8 time: 1.68
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
For this problem, if I set the number of Scenarios = 1, I am able to find a controller. It is when I try to append the additional LMI constraints and construct this controller using all random plant data (previously used for finding my R,S pair) that I encounter this dual infeasibility.
This boils down to my main questions:
1.) Does the infeasibility report indicate that I can bound this SDP to stop at a primal feasible solution? From studying dual infeasibility, it seems that I may be encountering a situation like
min x
s.t. x <= 5
Meaning my problem is unbounded, but I am running through feasible solutions toward -inf.
2.) If 1.) is the case, do you have a suggestion for bounding this problem using constructs or commands within YALMIP or MOSEK?
3.) Is the infeasibility so strong that my feasible space is actually null? Which output, here, would indicate this?
Thank you,
Chris
Johan Löfberg
Apr 11, 2018, 10:00:47 AM4/11/18
to YALMIP
If you want to understand if an infeasibility report really means infeasibility in your world (as it could relate to a dualized version etc after various reformulations), you simply remove the objective. If feasible, it is really an issue about unboundedness
These second recovery steps can be sensitive as they might assume that the first step was solved strictly. This is not guaranteed, as solvers typically generate slightly infeasible solutions. To combat that, you would have to add some margins in the initial step to ensure that the data you compute really is feasible
These second recovery steps can be sensitive as they might assume that the first step was solved strictly. This is not guaranteed, as solvers typically generate slightly infeasible solutions. To combat that, you would have to add some margins in the initial step to ensure that the data you compute really is feasible
Chris D'Angelo
Apr 11, 2018, 10:15:17 AM4/11/18
to YALMIP
I see. My first step involves solving feasibility problem, as well, as I search for the existence of R,S that satisfy LMI constraints. The code above also solves the feasibility problem - so I take it that you're saying that the issue does not appear to be unboundedness.
That leads me to your second point. What do you mean by "margins"? How would this be approached?
My feasible R, S pairs are not very close to being singular. Indeed, their eigenvalues are:
>> [eig(R) eig(S)]
ans =
8.8380 3.2631
9.5505 5.9450
9.6814 6.2135
9.7715 7.5082
9.8226 9.2184
9.8685 9.5792
9.8988 9.5988
9.9220 9.7131
9.9363 9.9193
9.9451 9.9506
9.9667 9.9556
9.9806 9.9628
10.0398 9.9634
10.0945 9.9687
10.2339 9.9722
10.3100 9.9778
10.3637 9.9832
10.9557 10.0292
12.2670 10.7684
12.2704 11.3171
Is there another indicator that I can use that gauges the quality of my solution? For only 10 scenarios, here is my solver output:
>> ScenarioControllerConstruction
Problem
Name :
Objective sense : min
Type : CONIC (conic optimization problem)
Constraints : 420
Cones : 0
Scalar variables : 0
Matrix variables : 21
Integer variables : 0
Optimizer started.
Presolve started.
Linear dependency checker started.
Linear dependency checker terminated.
Eliminator - tries : 0 time : 0.00
Lin. dep. - tries : 1 time : 0.00
Lin. dep. - number : 0
Presolve terminated. Time: 0.00
Problem
Name :
Objective sense : min
Type : CONIC (conic optimization problem)
Constraints : 420
Cones : 0
Scalar variables : 0
Matrix variables : 21
Integer variables : 0
Optimizer - threads : 8
Optimizer - solved problem : the primal
Optimizer - Constraints : 420
Optimizer - Cones : 0
Optimizer - Scalar variables : 0 conic : 0
Optimizer - Semi-definite variables: 21 scalarized : 6130
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 : 8.84e+04 after factor : 8.84e+04
Factor - dense dim. : 0 flops : 2.41e+08
ITE PFEAS DFEAS GFEAS PRSTATUS POBJ DOBJ MU TIME
0 1.8e+01 4.0e+00 5.0e+02 0.00e+00 1.998765864e+03 0.000000000e+00 1.0e+00 0.08
1 4.5e+00 9.7e-01 9.1e+01 -6.01e-01 1.229609824e+03 0.000000000e+00 2.4e-01 0.22
2 1.1e-01 2.4e-02 3.3e+00 1.73e-02 5.645740970e+01 0.000000000e+00 6.1e-03 0.34
3 7.3e-07 1.6e-07 7.5e-03 9.55e-01 3.699567407e-04 0.000000000e+00 4.0e-08 0.45
4 4.5e-14 2.8e-12 1.2e-12 1.00e+00 1.831652553e-11 0.000000000e+00 2.4e-15 0.54
Optimizer terminated. Time: 0.57
Interior-point solution summary
Problem status : PRIMAL_AND_DUAL_FEASIBLE
Solution status : OPTIMAL
Primal. obj: 1.8316525526e-11 nrm: 6e-13 Viol. con: 4e-10 barvar: 5e-22
Dual. obj: 0.0000000000e+00 nrm: 2e+01 Viol. con: 0e+00 barvar: 4e-11
Optimizer summary
Optimizer - time: 0.57
Interior-point - iterations : 4 time: 0.55
Chris D'Angelo
Apr 11, 2018, 10:34:47 AM4/11/18
to YALMIP
I now believe that I see what you mean.
I even went and added margin to R and S s.t. they are >=epsilon*eye(n), where epsilon is a small "tunable" number. I still found feasible R,S pairs with eigenvalues "far enough away from zero".
I then retooled my controller construction LMI to be:
Fcontroller = 0.5*(Fcontroller+Fcontroller')<=-.1*eye(mm);
where mm is the size of this LMI. I still encountered a dual infeasibility certificate. I'm trying to rule out numerical issues and appreciate your help with modeling tips and approaches. Can you see anything that I am missing, here, or can you suggest somewhere else for me to look?
Johan Löfberg
Apr 11, 2018, 11:23:47 AM4/11/18
to YALMIP
without a minimal reproducible example it is impossible to say anything
Johan Löfberg
Apr 11, 2018, 11:29:23 AM4/11/18
to YALMIP
psd:ness of R and S is not the hard part I guess. I presume you have [R I;I S] >=0, so that is the constraint that likely isn't definite enough (just a guess of course)
Johan Löfberg
Apr 11, 2018, 2:07:57 PM4/11/18
to YALMIP
Is the magic constant 0.1 known to be valid?
0.5*(Fcontroller+Fcontroller')<=-.1*eye(mm)
With the magic value 0.1, already optimize(G(3)) is infeasible. Changing the magic value 0.1 to 0 leads to all individual constraint being feasible, although already optimize(G(1:2)) is infeasible, so working with the case of 10 systems is not necessary in your debugging
btw, to avoid numerical issues causing lack of symmetry which forces you to symmetrize it, you should compute PSI_Xcl as
PSI_Xcl = [(Xcl*A0)'+Xcl*A0 Xcl*B0 C0';...
(Xcl*B0)' -gamma*eye(size(Xcl*B0,2)) D11';...
C0 D11 -gamma*eye(size(D11,1))];the first problem you solve appears very simple and the solution is far from feasibility border, hence non-strict solution there is definitely not the issue
My guess is a bug in your code (or theory), not really anything solver/yalmip related
Chris D'Angelo
Apr 11, 2018, 2:24:03 PM4/11/18
to YALMIP
The magic constant was arbitrarily chosen since I wanted to force this inequality to be strict (per our discussions above).
Thank you for the tip about calculation of PSI_Xcl.
I was thinking about this problem a bit more during a walk at lunch. I think that your guess about there being a problem with my theory is correct. Controller construction is dependent upon the plant data - finding feasible R,S pairs indicate the existence of controllers that meet one specified performance requirement for each scenario in the first composite LMI. So...I have basically "invented" (although it has probably been done before) a way to gain-schedule H-infinity controllers that produce the same closed-loop performance within a random plant set. Hence, our ability to synthesize a controller around the individual plant data.
Thank you very much for your help. I value your opinions and expertise (I've read several of your papers, as well).
Here's to research!
Reply all
Reply to author
Forward
0 new messages