|
ILMI algorithm 10 Months, 4 Weeks ago
|
Karma: 0
|
|
Iterative linear matrix inequality (ILMI) algorithm
Given & (with realization ( A, B, C)) stabilizable via state feed-back.
Step 1. Select Q>0, and solve P from the following algebraic
Riccati equation
A'P+PA-PBB'P+Q=0.
Set i=1 and X1=P.
Step 2. Solve the following optimization problem for Pi , F and ai.
OP1: Minimize ai subject to the following LMI constraints
[A'Pi+ Pi A-XiBB'Pi-PiBB'Xi+XiBB'Xi-aiPi (B'Pi+FC)' ; (B'Pi+FC) -I] < 0 (8)
Pi = Pi' > 0 (9)
Denote a*i as the minimized value of ai.
Step 3. If a*i <= 0, F is a stabilizing static output feedback gain. Stop.
Step 4. Solve the following optimization proble for Pi and F.
OP2: Minimize trace(Pi) subject to the above LMI constraints (8) and (9) with ai = a*i . Denote P*i as the Pi that mini-mized trace (Pi).
Step 5. If || Xi - P*i || < delta , a prescribed tolerance, go to Step 6,
else set i=i+1 and Xi = P*i-1 , then go to Step 2.
Step 6. The system may not be stabilizable via static output feedback. Stop
example :
A=[0 1;1 0] ; B=[1 0]' ; C=[1 15]
a=-0.0377 after 15 iterations and P=[9.2557 6.4888 ; 6.4848 92.9376] ; F=-0.7369
Could you please help me solving this problem in YALMIP ? Thank you!
|
|
|
|
|
|
|
Re: ILMI algorithm 10 Months, 4 Weeks ago
|
Karma: 32
|
|
Well, it is easily implemented, but how do intend to solve (8). It is bilinear in ai and Pi. Sure, easily solved using a bisection since it is quasi-convex, but are you aware of this. Your notation iterative LMI is not correct, it should be iterative bilinear matrix inequality, or iterative GEVP.
|
|
lofberg
Platinum Boarder
Posts: 2280
|
|
|
|
|
Re: ILMI algorithm 10 Months, 4 Weeks ago
|
Karma: 32
|
Solution (implements a generic GEVP solver)
| Code: |
function [P,F] = ilmi
Q = eye(2);
A=[0 1;1 0];
B=[1 0]';
C=[1 15];
delta = 1e-3;
sdpvar a
F = sdpvar(1,1);
P = sdpvar(2);
X = sdpvar(2);
X = are(A,B*B',Q);
i = 1;
while i <= 100
RHS = [P zeros(2,1); zeros(1,2) 0];
LHS = [A'*P+P*A-X*B*B'*P-P*B*B'*X+X*B*B'*X (B'*P+F*C)' ; (B'*P+F*C) -eye(1)];
% Find optimal ai
ai = solvelocalgevp(LHS,RHS,[],a)
% Now solve with that ai to recover the last feasible solution in the
% gevp code (could have been returned from there but it is easier to
% solve again...)
Constraints = [LHS <= ai*RHS,P >= 0];
solvesdp(Constraints,trace(P))
if ai <= 0
P = double(P);
F = double(F);
break
end
Pi = double(P);
if norm(X-Pi) <= delta
P = [];
F = [];
break
end
X = Pi;
end
function t_upper = solvelocalgevp(A,B,C,t)
% Solves min t s.t A(x)<=t*B(x), B>=0, C>=0
% Find initial lower and upper bounds
ops = sdpsettings('verbose',0);
sol = solvesdp([C>=0, A <= 0*B, B>=0],[],ops);
if sol.problem == 0
% Ok, 0 is an upper_bound, i.e. feasible
t_upper = 0;
t_lower = -1;
while sol.problem == 0
t_lower = (t_lower)*2;
sol = solvesdp([C>=0, A <= t_lower*B, B>=0],[],ops);
end
elseif sol.problem == 1
% Ok 0 is a lower bound, i.e. infeasible
t_lower = 0;
t_upper = 1;
while sol.problem == 1
t_upper = t_upper*2;
sol = solvesdp([C>=0, A <= t_upper*B, B>=0],[],ops);
end
else
error('Hmm, you''d better robustify the code!')
end
% OK, we have the bounds, start bisecting
while t_upper-t_lower >= 1e-3
t_mid = (t_upper+t_lower)/2;
sol = solvesdp([C>=0, A <= t_mid*B, B>=0],[],ops);
if sol.problem == 0
t_upper = t_mid;
else
t_lower = t_mid;
end
disp([t_lower t_upper]);
end
|
|
|
lofberg
Platinum Boarder
Posts: 2280
|
|
Last Edit: 2012/06/28 03:21 By lofberg.
|
|
|
Re: ILMI algorithm 10 Months, 4 Weeks ago
|
Karma: 0
|
|
Thank you Lofberg for the advices . In fact , it's SOF (Static Output Feedback) resolved by an iterative LMI (iterative GEVP ). But when I simulated the program and it showed , failure of FEASP . Perhaps I should change the parameters of the system .
|
|
|
|
|
|
|
|
lofberg
Platinum Boarder
Posts: 2280
|
|
|
|
|
Re: ILMI algorithm 10 Months, 4 Weeks ago
|
Karma: 0
|
|
I have obtained the following results by implementing the previous program :
SeDuMi 1.1 by AdvOL, 2005 and Jos F. Sturm, 1998, 2001-2003.
Alg = 2: xz-corrector, theta = 0.250, beta = 0.500
eqs m = 4, order n = 6, dim = 14, blocks = 3
nnz(A) = 14 + 0, nnz(ADA) = 16, nnz(L) = 10
it : b*y gap delta rate t/tP* t/tD* feas cg cg prec
0 : 5.62E+000 0.000
1 : -2.45E+000 5.12E-001 0.000 0.0911 0.9900 0.9900 1.38 1 1 9.1E-001
2 : -2.46E+000 1.11E-001 0.000 0.2170 0.9000 0.9000 1.21 1 1 1.8E-001
3 : -2.56E+000 2.99E-003 0.000 0.0270 0.9900 0.9900 0.94 1 1 5.0E-003
4 : -2.56E+000 5.30E-004 0.000 0.1768 0.9000 0.9000 1.01 1 1 8.8E-004
5 : -2.56E+000 3.20E-005 0.347 0.0603 0.9900 0.9900 1.01 1 1 5.2E-005
6 : -2.56E+000 7.82E-006 0.000 0.2447 0.9000 0.9000 1.00 1 1 1.3E-005
7 : -2.56E+000 2.49E-007 0.000 0.0319 0.9900 0.9900 1.00 1 1 4.0E-007
8 : -2.56E+000 4.54E-008 0.000 0.1821 0.9000 0.9000 1.00 2 2 7.3E-008
9 : -2.56E+000 2.61E-009 0.313 0.0575 0.9900 0.9900 1.00 2 2 4.2E-009
10 : -2.56E+000 1.19E-010 0.000 0.0455 0.9900 0.9900 1.00 2 2 1.9E-010
iter seconds digits c*x b*y
10 0.4 Inf -2.5582458114e+000 -2.5582458113e+000
|Ax-b| = 1.7e-010, [Ay-c]_+ = 4.9E-011, |x|= 2.7e+000, |y|= 2.0e+000
Detailed timing (sec)
Pre IPM Post
3.120E-002 4.212E-001 0.000E+000
Max-norms: ||b||=1, ||c|| = 2.488584e+000,
Cholesky |add|=0, |skip| = 0, ||L.L|| = 1.27425.
ans =
yalmiptime: 0.221000000000000
solvertime: 0.251000000000001
info: 'Successfully solved (SeDuMi-1.1)'
problem: 0
dimacs: [NaN NaN NaN NaN NaN NaN]
ans =
[]
How can I retrieve the values of P , F and ai ?
|
|
|
|
|
|
|
