Teljes Matlab script kiegészítő függvényekkel.
File: d2018_01_31_K_prelim_L_codesign_v6_rhofuggoBL_min.m Directory: projects/3_outsel/2017_11_13_lpv_passivity Author: Peter Polcz (ppolcz@gmail.com)
Created on 2018. January 31.
global VERBOSE SCOPE_DEPTH
VERBOSE = 1;
SCOPE_DEPTH = -1;
s = tf('s');
H = @(s) [
(s-1)/(s-2)/(s+1) 1/(s+3)/(s-0.1)
(s-7)/(s+1)/(s+5) (s-6)/(s^2+5*s+6)
];
sys = minreal( ss( H(s) ) );
[A0,B0,C,D] = deal(sys.a, sys.b, sys.c, sys.d);
tol = 1e-10;
prec = -log10(tol);
A0 = round(A0,prec);
B0 = round(B0,prec);
C = round(C,prec);
D = round(D,prec);
[POLES,ZEROS] = pzmap(sys)
2 states removed. POLES = -5.0000 -1.0000 2.0000 -3.0000 -2.0000 0.1000 ZEROS = 5.8894 + 0.0000i -4.3164 + 0.0000i 0.7635 + 0.7978i 0.7635 - 0.7978i
rho_lim = [
-1 1
];
rho_rand = @(varargin) rand(varargin{:})*(rho_lim(2) - rho_lim(1)) + rho_lim(1);
A1 = A0;
A1(abs(A0) < 1) = 0;
A1 = A1 .* randn(size(A1))/10;
B1 = [ B0(:,2)/2 B0(:,2)*0 ];
% _fh: `funcion handle`
A_fh = @(rho) A0 + rho*A1;
B_fh = @(rho) B0 + rho*B1;
% Dimenziok
n_x = size(A0,1)
n_u = size(B0,2)
n_y = size(C,1)
n_r = n_u
n_yp = n_r
n_x = 6 n_u = 2 n_y = 2 n_r = 2 n_yp = 2
Kvadratikus stabilitas
Q = sdpvar(n_x);
N = sdpvar(n_u,n_x,'full');
Big_M = @(rho) A_fh(rho)'*Q + Q*A_fh(rho) - B_fh(rho)*N - N'*B_fh(rho)';
I = eye(n_x);
Constraints = [
Q - I >= 0
Big_M(rho_lim(1)) + 10*I <= 0
Big_M(rho_lim(2)) + 10*I <= 0
];
optimize(Constraints)
N = value(N);
Q = value(Q);
K = N/Q;
G = eye(n_u);
% for rho = rho_lim(1):0.25:rho_lim(2)
% eig(A_fh(rho) - B_fh(rho)*K)
% end
Problem Name : Objective sense : min Type : CONIC (conic optimization problem) Constraints : 33 Cones : 0 Scalar variables : 0 Matrix variables : 3 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 : 33 Cones : 0 Scalar variables : 0 Matrix variables : 3 Integer variables : 0 Optimizer - threads : 4 Optimizer - solved problem : the primal Optimizer - Constraints : 33 Optimizer - Cones : 0 Optimizer - Scalar variables : 0 conic : 0 Optimizer - Semi-definite variables: 3 scalarized : 63 Factor - setup time : 0.00 dense det. time : 0.00 Factor - ML order time : 0.00 GP order time : 0.00 Factor - nonzeros before factor : 561 after factor : 561 Factor - dense dim. : 0 flops : 5.11e+04 ITE PFEAS DFEAS GFEAS PRSTATUS POBJ DOBJ MU TIME 0 2.4e+00 6.0e+00 6.2e+01 0.00e+00 -2.520000000e+02 0.000000000e+00 1.0e+00 0.00 1 3.5e-01 9.0e-01 1.8e+00 -1.28e+00 -8.483253758e+02 0.000000000e+00 1.5e-01 0.01 2 5.9e-02 1.5e-01 3.6e-01 -1.44e+00 -2.195417908e+02 0.000000000e+00 2.5e-02 0.01 3 5.8e-04 1.5e-03 6.0e-02 8.47e-01 -1.774232594e+00 0.000000000e+00 2.4e-04 0.01 4 1.5e-12 3.7e-12 3.8e-11 9.96e-01 -6.194005287e-09 0.000000000e+00 6.1e-13 0.01 Optimizer terminated. Time: 0.01 Interior-point solution summary Problem status : PRIMAL_AND_DUAL_FEASIBLE Solution status : OPTIMAL Primal. obj: -6.1940052866e-09 nrm: 2e-10 Viol. con: 7e-10 barvar: 1e-18 Dual. obj: 0.0000000000e+00 nrm: 1e+02 Viol. con: 0e+00 barvar: 2e-10 Optimizer summary Optimizer - time: 0.01 Interior-point - iterations : 4 time: 0.01 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 ans = struct with fields: yalmiptime: 0.8212 solvertime: 0.0209 info: 'Successfully solved (MOSEK)' problem: 0
C1 = sdpvar([n_yp n_yp], [n_x n_x],'full');
D1 = sdpvar([n_yp n_yp], [n_y n_y],'full');
Q = sdpvar(n_x,n_x,'symmetric');
S = sdpvar(n_x,n_x,'symmetric');
N = sdpvar([n_x n_x],[n_y n_y],'full');
P = blkdiag(Q,S);
C1_fh = @(rho) C1{1} + rho*C1{2};
D1_fh = @(rho) D1{1} + rho*D1{2};
N_fh = @(rho) N{1} + rho*N{2};
% Open loop matrices
Bo_fh = @(rho) [
B_fh(rho)
zeros(n_x,n_u)
];
Co_fh = @(rho) [
D1_fh(rho)*C+C1_fh(rho) -C1_fh(rho)
];
% Closed loop matrices
Bc_fh = @(rho) Bo_fh(rho)*G;
AcP_PAc_fh = @(rho) [
Q*(A_fh(rho)-B_fh(rho)*K) + (A_fh(rho)-B_fh(rho)*K)'*Q , Q*B_fh(rho)*K
K'*B_fh(rho)'*Q , S*A_fh(rho) + A_fh(rho)'*S - N_fh(rho)*C - C'*N_fh(rho)'
];
W = eye(n_yp);
Lambda2 = @(rho) [
AcP_PAc_fh(rho) , P*Bc_fh(rho)-Co_fh(rho)' , Co_fh(rho)'
Bc_fh(rho)'*P-Co_fh(rho) , zeros(n_r,n_r) , zeros(n_r,n_yp)
Co_fh(rho) , zeros(n_yp,n_r) , -inv(W)
];
Constraints = [
Lambda2(rho_lim(1)) <= 0
Lambda2(rho_lim(2)) <= 0
P - 0.0001*eye(size(P)) >= 0
]
++++++++++++++++++++++++++++++++++ | ID| Constraint| ++++++++++++++++++++++++++++++++++ | #1| Matrix inequality 16x16| | #2| Matrix inequality 16x16| | #3| Matrix inequality 12x12| ++++++++++++++++++++++++++++++++++
sdpopts = sdpsettings('solver','sedumi');
optimize(Constraints,[],sdpopts)
check(Constraints)
The coefficient matrix is not full row rank, numerical problems may occur. SeDuMi 1.3 by AdvOL, 2005-2008 and Jos F. Sturm, 1998-2003. Alg = 2: xz-corrector, theta = 0.250, beta = 0.500 eqs m = 98, order n = 45, dim = 657, blocks = 4 nnz(A) = 1962 + 0, nnz(ADA) = 9604, nnz(L) = 4851 it : b*y gap delta rate t/tP* t/tD* feas cg cg prec 0 : 3.12E+01 0.000 1 : 0.00E+00 8.38E+00 0.000 0.2687 0.9000 0.9000 2.00 1 1 9.0E+00 2 : 0.00E+00 2.97E+00 0.000 0.3548 0.9000 0.9000 1.27 1 1 3.3E+00 3 : 0.00E+00 1.27E+00 0.000 0.4257 0.9000 0.9000 1.09 1 1 1.7E+00 4 : 0.00E+00 6.36E-01 0.000 0.5023 0.9000 0.9000 1.04 1 1 1.1E+00 5 : 0.00E+00 2.05E-01 0.000 0.3221 0.9000 0.9000 1.01 1 1 7.6E-01 6 : 0.00E+00 8.58E-02 0.000 0.4190 0.9000 0.9000 1.00 1 1 6.6E-01 7 : 0.00E+00 4.02E-02 0.000 0.4680 0.9000 0.9000 0.98 1 1 6.3E-01 8 : 0.00E+00 1.46E-02 0.000 0.3646 0.9000 0.9000 0.98 1 1 6.1E-01 9 : 0.00E+00 5.98E-03 0.000 0.4084 0.9000 0.9000 0.97 1 1 3.3E-01 10 : 0.00E+00 1.97E-03 0.000 0.3294 0.9000 0.9000 0.96 1 1 1.1E-01 11 : 0.00E+00 7.36E-04 0.000 0.3737 0.9000 0.9000 0.94 1 1 4.6E-02 12 : 0.00E+00 2.51E-04 0.000 0.3406 0.9000 0.9000 0.91 1 1 1.7E-02 13 : 0.00E+00 9.50E-05 0.000 0.3788 0.9000 0.9000 0.84 1 1 7.5E-03 14 : 0.00E+00 3.38E-05 0.000 0.3555 0.9000 0.9000 0.75 1 1 3.3E-03 15 : 0.00E+00 1.36E-05 0.000 0.4019 0.9000 0.9000 0.59 1 1 1.9E-03 16 : 0.00E+00 4.78E-06 0.000 0.3519 0.9000 0.9000 0.44 1 1 9.3E-04 17 : 0.00E+00 1.90E-06 0.000 0.3972 0.9000 0.9000 0.27 1 1 3.9E-04 18 : 0.00E+00 6.49E-07 0.000 0.3419 0.9000 0.9000 0.18 1 1 3.3E-06 19 : 0.00E+00 2.48E-07 0.000 0.3826 0.9000 0.9000 0.09 1 1 2.0E-06 20 : 0.00E+00 8.49E-08 0.000 0.3421 0.9000 0.9000 0.07 1 1 1.1E-06 21 : 0.00E+00 3.23E-08 0.000 0.3802 0.9000 0.9000 0.02 1 1 7.0E-07 22 : 0.00E+00 1.11E-08 0.000 0.3446 0.9000 0.9000 0.04 1 1 3.9E-07 23 : 0.00E+00 4.26E-09 0.000 0.3831 0.9000 0.9000 -0.01 1 1 2.5E-07 24 : 0.00E+00 1.47E-09 0.000 0.3447 0.9000 0.9000 0.02 1 1 1.4E-07 25 : 0.00E+00 5.63E-10 0.000 0.3835 0.9000 0.9000 -0.01 1 1 9.2E-08 26 : 0.00E+00 1.94E-10 0.000 0.3437 0.9000 0.9000 0.02 1 1 5.1E-08 27 : 0.00E+00 7.40E-11 0.000 0.3823 0.9000 0.9000 -0.02 1 2 3.3E-08 28 : 0.00E+00 2.54E-11 0.000 0.3435 0.9000 0.9000 0.01 1 1 1.9E-08 29 : 0.00E+00 9.71E-12 0.000 0.3818 0.9000 0.9000 -0.02 1 2 1.2E-08 30 : 0.00E+00 3.33E-12 0.000 0.3436 0.9000 0.9000 0.01 2 2 6.7E-09 31 : 0.00E+00 1.27E-12 0.000 0.3818 0.9000 0.9000 -0.02 2 2 4.4E-09 32 : 0.00E+00 4.37E-13 0.000 0.3434 0.9000 0.9000 0.01 2 2 2.4E-09 33 : 0.00E+00 1.67E-13 0.000 0.3817 0.9000 0.9000 -0.02 2 2 1.6E-09 34 : 0.00E+00 5.73E-14 0.000 0.3434 0.9000 0.9000 0.01 2 2 8.8E-10 iter seconds digits c*x b*y 34 0.4 11.4 -4.0581649749e-07 0.0000000000e+00 |Ax-b| = 8.8e-10, [Ay-c]_+ = 3.8E-11, |x|= 9.7e+04, |y|= 3.7e+05 Detailed timing (sec) Pre IPM Post 8.228E-03 1.423E-01 2.303E-03 Max-norms: ||b||=0, ||c|| = 1, Cholesky |add|=0, |skip| = 2, ||L.L|| = 8.17022. ans = struct with fields: yalmiptime: 0.1423 solvertime: 0.1536 info: 'Numerical problems (SeDuMi-1.3)' problem: 4 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ | ID| Constraint| Primal residual| Dual residual| ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ | #1| Matrix inequality| -2.9371e-11| 9.0952e-12| | #2| Matrix inequality| -3.888e-11| 9.0296e-12| | #3| Matrix inequality| 3.883e-05| 1.3052e-12| ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Q = value(Q)
S = value(S)
P = blkdiag(Q,S);
N = cellfun(@value, N, 'UniformOutput', false);
C1 = cellfun(@value, C1, 'UniformOutput', false);
D1 = cellfun(@value, D1, 'UniformOutput', false);
L = cellfun(@(N) S\N, N, 'UniformOutput', false);
L_fh = @(rho) L{1} + rho*L{2};
C1_fh = @(rho) C1{1} + rho*C1{2};
D1_fh = @(rho) D1{1} + rho*D1{2};
Q = 1.0e-03 * 0.2826 0.0329 -0.1081 0.0319 0.0520 0.0758 0.0329 0.3330 -0.1131 0.0654 0.1099 0.0341 -0.1081 -0.1131 0.2819 0.0131 -0.1024 -0.0241 0.0319 0.0654 0.0131 0.5366 0.0444 0.0160 0.0520 0.1099 -0.1024 0.0444 0.3027 -0.0164 0.0758 0.0341 -0.0241 0.0160 -0.0164 0.3959 S = 1.0e+05 * 0.1291 0.0768 -0.0130 0.1568 0.0181 -0.1400 0.0768 0.2639 0.0221 -0.0019 -0.0402 -0.0478 -0.0130 0.0221 0.1189 0.2056 -0.0220 0.0366 0.1568 -0.0019 0.2056 2.4612 -0.0897 -0.7398 0.0181 -0.0402 -0.0220 -0.0897 0.0680 -0.0376 -0.1400 -0.0478 0.0366 -0.7398 -0.0376 0.4084\begin{align} {\LARGE(12) \quad} L(\rho) = \left(\begin{array}{cc} 3.92\rho -32.8 & 0.676\rho -10.7 \\ 6.76\rho -47.7 & 1.37\rho -18.6 \\ 120.0-17.0\rho & 46.3-3.51\rho \\ 7.8\rho -54.3 & 1.53\rho -21.0 \\ 20.0\rho -148.0 & 3.95\rho -56.5 \\ 19.6\rho -137.0 & 3.84\rho -53.3 \\ \end{array}\right) \end{align} \begin{align} {\LARGE(13) \quad} C_1(\rho) = \left(\begin{array}{cccccc} 2.0310^{-4}-8.7510^{-5}\rho & 1.810^{-5}\rho -1.5710^{-4} & 1.8210^{-4}\rho +2.7110^{-4} & 9.1110^{-6}\rho -2.0410^{-4} & 1.8110^{-4}\rho -4.2310^{-4} & 6.7610^{-6}\rho -1.8910^{-4} \\ 2.5310^{-9}\rho -1.7510^{-4} & 1.3610^{-8}\rho +3.6110^{-5} & 7.510^{-9}\rho +3.6410^{-4} & 2.3610^{-8}\rho +1.8210^{-5} & 3.6210^{-4}-8.7110^{-9}\rho & 1.3510^{-5}-3.8310^{-9}\rho \\ \end{array}\right) \end{align} \begin{align} {\LARGE(14) \quad} D_1(\rho) = \left(\begin{array}{cc} 4.7810^{-5}\rho +8.6310^{-5} & 3.3210^{-5}\rho +9.8710^{-5} \\ 9.5610^{-5}-2.5110^{-8}\rho & 2.4410^{-9}\rho +6.6410^{-5} \\ \end{array}\right) \end{align}
Ao_fh = @(rho) [
A_fh(rho) zeros(n_x,n_x)
zeros(n_x,n_x) A_fh(rho)-L_fh(rho)*C
];
Bo_fh = @(rho) [
B_fh(rho)
zeros(n_x,n_u)
];
Co_fh = @(rho) [
D1_fh(rho)*C+C1_fh(rho) -C1_fh(rho)
];
Do = zeros(n_yp,n_r);
% Closed loop matrices
Ac_fh = @(rho) [
A_fh(rho)-B_fh(rho)*K B_fh(rho)*K
zeros(n_x,n_x) A_fh(rho)-L_fh(rho)*C
];
Bc_fh = @(rho) Bo_fh(rho)*G;
for rho = rho_lim(1):0.25:rho_lim(2)
sys = ss(Ao_fh(rho),Bo_fh(rho),Co_fh(rho),Do);
[POLES,ZEROS] = pzmap(sys);
fprintf('\nzeros if rho = %g\n', rho)
disp(ZEROS)
pcz_info(all(real(ZEROS) < 0), 'Zeros are all negative');
end
zeros if rho = -1 -13.7554 + 9.9143i -13.7554 - 9.9143i -4.0327 + 0.0000i -1.2935 + 1.5082i -1.2935 - 1.5082i -0.4362 + 0.0000i -1.1576 + 0.9804i -1.1576 - 0.9804i -0.2595 + 0.0000i -0.6534 + 0.0000i [ OK ] Zeros are all negative zeros if rho = -0.75 -13.3061 + 9.8889i -13.3061 - 9.8889i -3.9972 + 0.0000i -1.3292 + 1.4350i -1.3292 - 1.4350i -0.4480 + 0.0000i -1.1472 + 0.9825i -1.1472 - 0.9825i -0.2945 + 0.0000i -0.5798 + 0.0000i [ OK ] Zeros are all negative zeros if rho = -0.5 -12.8542 + 9.8598i -12.8542 - 9.8598i -3.9621 + 0.0000i -1.3654 + 1.3495i -1.3654 - 1.3495i -1.1367 + 0.9848i -1.1367 - 0.9848i -0.4639 + 0.0000i -0.3647 + 0.0000i -0.4711 + 0.0000i [ OK ] Zeros are all negative zeros if rho = -0.25 -12.3992 + 9.8273i -12.3992 - 9.8273i -3.9276 + 0.0000i -1.4013 + 1.2478i -1.4013 - 1.2478i -0.4856 + 0.0000i -1.1261 + 0.9872i -1.1261 - 0.9872i -0.3988 + 0.1174i -0.3988 - 0.1174i [ OK ] Zeros are all negative zeros if rho = 0 -11.9410 + 9.7915i -11.9410 - 9.7915i -3.8939 + 0.0000i -1.4358 + 1.1232i -1.4358 - 1.1232i -0.5159 + 0.0000i -1.1153 + 0.9896i -1.1153 - 0.9896i -0.3799 + 0.1716i -0.3799 - 0.1716i [ OK ] Zeros are all negative zeros if rho = 0.25 -11.4791 + 9.7526i -11.4791 - 9.7526i -3.8614 + 0.0000i -1.4662 + 0.9632i -1.4662 - 0.9632i -0.5605 + 0.0000i -1.1044 + 0.9921i -1.1044 - 0.9921i -0.3611 + 0.2104i -0.3611 - 0.2104i [ OK ] Zeros are all negative zeros if rho = 0.5 -11.0131 + 9.7111i -11.0131 - 9.7111i -3.8307 + 0.0000i -1.4852 + 0.7353i -1.4852 - 0.7353i -0.6343 + 0.0000i -1.0933 + 0.9946i -1.0933 - 0.9946i -0.3425 + 0.2412i -0.3425 - 0.2412i [ OK ] Zeros are all negative zeros if rho = 0.75 -10.5427 + 9.6673i -10.5427 - 9.6673i -3.8027 + 0.0000i -1.4491 + 0.2060i -1.4491 - 0.2060i -0.8246 + 0.0000i -1.0820 + 0.9972i -1.0820 - 0.9972i -0.3240 + 0.2668i -0.3240 - 0.2668i [ OK ] Zeros are all negative zeros if rho = 1 -10.0676 + 9.6218i -10.0676 - 9.6218i -3.7794 + 0.0000i -2.2157 + 0.0000i -0.8149 + 0.3839i -0.8149 - 0.3839i -1.0706 + 0.9997i -1.0706 - 0.9997i -0.3057 + 0.2888i -0.3057 - 0.2888i [ OK ] Zeros are all negative
concat = @(fh) { fh(0) , fh(1)-fh(0) };
A = concat(Ao_fh);
B = concat(Bo_fh);
C = concat(Co_fh);
D = {Do, Do*0};
AA = [A{:}];
BB = [B{:}];
Im_B = IMA(BB);
L = ORTCO(Im_B);
tol = 1e-5;
Ker_C = INTS(KER(C{1}), KER(C{2}));
Ker_C = pcz_vecalg_ints(pcz_vecalg_null(C{1},tol), pcz_vecalg_null(C{2},tol),tol);
[R,V] = CSA(AA, Im_B, Ker_C);
dim_V = rank(V);
dim_Ker_V = size(V,1) - dim_V;
dim_L = rank(L);
assert(rank(R) == 0 && rank([V Ker_C]) == rank(Ker_C), ...
'Strong invertibility condition is not satisfied!');
assert(rank(INTS(V, Im_B)) == 0)
assert(dim_L >= dim_V);
T = [ ORTCO(V) L(:,1:dim_V) ]';
N = numel(A);
tA = cell(1,N);
tB = cell(1,N);
tC = cell(1,N);
tD = cell(1,N);
for i = 1:numel(A)
tA{i} = T*A{i}/T;
tB{i} = T*B{i};
tC{i} = C{i}/T;
tD{i} = D{i};
end
tA_fh = @(rho) tA{1} + tA{2}*rho;
tB_fh = @(rho) tB{1} + tB{2}*rho;
tC_fh = @(rho) tC{1} + tC{2}*rho;
tD_fh = @(rho) tD{1} + tD{2}*rho;
\begin{align} {\LARGE(15) \quad}
\bar A(\rho) =
\left(\begin{array}{cccccccccccc|}
0.24\rho -5.4 & -0.93 & -0.46 & 0.093 & 0.46 & 0.46 & 0 & 0 & 0 & 0 & 0 & 0 \\
1.9-0.029\rho & -0.32 & -0.19 & 0.038 & 0.19 & 0.19 & 0 & 0 & 0 & 0 & 0 & 0 \\
-0.29 & -0.73 & -0.17\rho -1.4 & 0.29 & 4.810^{-3}\rho -1.5 & 0.35\rho -1.5 & 0 & 0 & 0 & 0 & 0 & 0 \\
0.059 & 0.15 & 0.22 & 0.055 & 0.28 & 0.28 & 0 & 0 & 0 & 0 & 0 & 0 \\
0.29 & 0.73 & -0.3\rho -3.2 & 0.64 & -0.077\rho -1.8 & 0.19 & 0 & 0 & 0 & 0 & 0 & 0 \\
0.29 & 0.73 & 0.1 & -0.02 & 1.9-0.071\rho & -0.1 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 7.7-0.9\rho & -0.81\rho -0.41 & 3.2-0.7\rho & 32.0-3.8\rho & 0.36\rho +2.2 & 1.7\rho -19.0 \\
0 & 0 & 0 & 0 & 0 & 0 & 23.0-2.1\rho & -1.1\rho -3.4 & 3.7-1.1\rho & 47.0-6.5\rho & 0.44\rho +5.6 & 3.2\rho -32.0 \\
0 & 0 & 0 & 0 & 0 & 0 & 5.3\rho -53.0 & 2.8\rho +6.4 & 2.6\rho -11.0 & 16.0\rho -122.0 & -1.0\rho -15.0 & 78.0-7.7\rho \\
0 & 0 & 0 & 0 & 0 & 0 & 24.0-2.4\rho & -1.4\rho -3.1 & 4.8-1.3\rho & 53.0-7.5\rho & 0.55\rho +6.2 & 3.6\rho -36.0 \\
0 & 0 & 0 & 0 & 0 & 0 & 65.0-6.1\rho & -3.5\rho -7.4 & 9.5-3.6\rho & 166.0-19.0\rho & 1.3\rho +14.0 & 9.3\rho -97.0 \\
0 & 0 & 0 & 0 & 0 & 0 & 61.0-6.0\rho & -3.5\rho -8.0 & 11.0-3.3\rho & 133.0-19.0\rho & 1.3\rho +17.0 & 9.1\rho -91.0 \\ \hline
\end{array}\right)
\end{align}
\begin{align} {\LARGE(16) \quad}
\bar B(\rho) =
\left(\begin{array}{cc}
2.0 & 0 \\
0 & 0 \\
\rho +1.1 & 2.0 \\
-0.23 & 0 \\
\rho -1.1 & 2.0 \\
-1.1 & 0 \\
0 & 0 \\
0 & 0 \\
0 & 0 \\
0 & 0 \\
0 & 0 \\
0 & 0 \\ \hline
\end{array}\right)
\end{align}
\begin{align} {\LARGE(17) \quad}
\bar C(\rho) =
\left(\begin{array}{cccccccccccc|}
2.910^{-4}-5.610^{-5}\rho & -3.210^{-6}\rho -2.410^{-4} & 1.810^{-4}\rho +2.510^{-4} & 5.710^{-5}\rho -1.110^{-4} & 2.010^{-4}\rho -3.510^{-4} & -4.110^{-5}\rho -3.110^{-4} & 8.810^{-5}\rho -2.010^{-4} & 1.610^{-4}-1.810^{-5}\rho & -1.810^{-4}\rho -2.710^{-4} & 2.010^{-4}-9.110^{-6}\rho & 4.210^{-4}-1.810^{-4}\rho & 1.910^{-4}-6.810^{-6}\rho \\
2.110^{-12}\rho -1.110^{-4} & 4.410^{-12}\rho -6.410^{-6} & 3.610^{-4}-4.710^{-12}\rho & 2.910^{-12}\rho +1.110^{-4} & 9.310^{-12}\rho +4.010^{-4} & 5.510^{-12}\rho -8.110^{-5} & 1.710^{-4}-2.510^{-9}\rho & -1.410^{-8}\rho -3.610^{-5} & -7.510^{-9}\rho -3.610^{-4} & -2.410^{-8}\rho -1.810^{-5} & 8.710^{-9}\rho -3.610^{-4} & 3.810^{-9}\rho -1.410^{-5} \\
\end{array}\right)
\end{align}
(Ennek most nincs jelentosege)
tP = T' \ P / T;
tP22 = tP(dim_Ker_V+1:end,dim_Ker_V+1:end);
for rho = rho_lim
tA_rho = tA_fh(rho);
tA22 = tA_rho(dim_Ker_V+1:end,dim_Ker_V+1:end);
eig(tP22 * tA22 + tA22' * tP22)
end
ans = 0×1 empty double column vector ans = 0×1 empty double column vector
Using the output zeroing input
syms rho drho
A_sym = Ao_fh(rho);
B_sym = Bo_fh(rho);
C_sym = Co_fh(rho);
D_sym = Do;
dC_sym = Co_fh(drho);
A_zd = A_sym - B_sym*((C_sym*B_sym)\(C_sym*A_sym + dC_sym));
A_zd_fh = matlabFunction(A_zd, 'vars', [rho,drho]);
IS_OK = 1;
for rho = rho_lim(1):0.25:rho_lim(2)
for drho = -10:5:10
if pcz_info(all(real(eig(A_zd_fh(rho,drho))) < 0), ...
'Zero dynamics is stable if rho = %g, drho = %g.', rho, drho)
A_zd_eigvals = eig(A_zd_fh(rho,drho))
IS_OK = 0;
end
end
end
pcz_info(IS_OK, 'The zero dynamics is stable in everry sampling point')
[ OK ] Zero dynamics is stable if rho = -1, drho = -10. [ OK ] Zero dynamics is stable if rho = -1, drho = -5. [ OK ] Zero dynamics is stable if rho = -1, drho = 0. [ OK ] Zero dynamics is stable if rho = -1, drho = 5. [ OK ] Zero dynamics is stable if rho = -1, drho = 10. [ OK ] Zero dynamics is stable if rho = -0.75, drho = -10. [ OK ] Zero dynamics is stable if rho = -0.75, drho = -5. [ OK ] Zero dynamics is stable if rho = -0.75, drho = 0. [ OK ] Zero dynamics is stable if rho = -0.75, drho = 5. [ OK ] Zero dynamics is stable if rho = -0.75, drho = 10. [ OK ] Zero dynamics is stable if rho = -0.5, drho = -10. [ OK ] Zero dynamics is stable if rho = -0.5, drho = -5. [ OK ] Zero dynamics is stable if rho = -0.5, drho = 0. [ OK ] Zero dynamics is stable if rho = -0.5, drho = 5. [ OK ] Zero dynamics is stable if rho = -0.5, drho = 10. [ OK ] Zero dynamics is stable if rho = -0.25, drho = -10. [ OK ] Zero dynamics is stable if rho = -0.25, drho = -5. [ OK ] Zero dynamics is stable if rho = -0.25, drho = 0. [ OK ] Zero dynamics is stable if rho = -0.25, drho = 5. [ OK ] Zero dynamics is stable if rho = -0.25, drho = 10. [ OK ] Zero dynamics is stable if rho = 0, drho = -10. [ OK ] Zero dynamics is stable if rho = 0, drho = -5. [ OK ] Zero dynamics is stable if rho = 0, drho = 0. [ OK ] Zero dynamics is stable if rho = 0, drho = 5. [ OK ] Zero dynamics is stable if rho = 0, drho = 10. [ OK ] Zero dynamics is stable if rho = 0.25, drho = -10. [ OK ] Zero dynamics is stable if rho = 0.25, drho = -5. [ OK ] Zero dynamics is stable if rho = 0.25, drho = 0. [ OK ] Zero dynamics is stable if rho = 0.25, drho = 5. [ OK ] Zero dynamics is stable if rho = 0.25, drho = 10. [ OK ] Zero dynamics is stable if rho = 0.5, drho = -10. [ OK ] Zero dynamics is stable if rho = 0.5, drho = -5. [ OK ] Zero dynamics is stable if rho = 0.5, drho = 0. [ OK ] Zero dynamics is stable if rho = 0.5, drho = 5. [ OK ] Zero dynamics is stable if rho = 0.5, drho = 10. [ OK ] Zero dynamics is stable if rho = 0.75, drho = -10. [ OK ] Zero dynamics is stable if rho = 0.75, drho = -5. [ OK ] Zero dynamics is stable if rho = 0.75, drho = 0. [ OK ] Zero dynamics is stable if rho = 0.75, drho = 5. [ OK ] Zero dynamics is stable if rho = 0.75, drho = 10. [ OK ] Zero dynamics is stable if rho = 1, drho = -10. [ OK ] Zero dynamics is stable if rho = 1, drho = -5. [ OK ] Zero dynamics is stable if rho = 1, drho = 0. [ OK ] Zero dynamics is stable if rho = 1, drho = 5. [ OK ] Zero dynamics is stable if rho = 1, drho = 10. [ OK ] The zero dynamics is stable in everry sampling point
tol = 1e-10;
prec = -log10(tol);
is_negdef = @(A) all( round(eig(A), prec) <= 0 )
is_negdef( Ao_fh(rho_lim(1))'*P + P*Ao_fh(rho_lim(1)) )
is_negdef( Ao_fh(rho_lim(2))'*P + P*Ao_fh(rho_lim(2)) )
is_negdef = function_handle with value: @(A)all(round(eig(A),prec)<=0) ans = logical 0 ans = logical 0