Teljes Matlab script kiegészítő függvényekkel.
File: d2018_02_07_K_prelim_L_codesign_v6_stabilitas_hataran.m Directory: 1_projects/3_outsel/2017_11_13_lpv_passivity Author: Peter Polcz (ppolcz@gmail.com)
Created on 2018. February 07.
Inherited from
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 SCOPE_DEPTH VERBOSE LATEX_EQNR
SCOPE_DEPTH = 0;
VERBOSE = 1;
LATEX_EQNR = 0;
s = tf('s');
% NEM JO
H = @(s) [
(s-2) -3*(s+1)
s s^2-1
] / ( (s+2)*(s^2+4)*(s^2-s+1));
% JO
H = @(s) (s^2+1)/(s^2+3*s+2)/(s+5);
% EZ IS JO
H = @(s) [
s-2 -3
s s-1
] / (s^2+3*s+2)/(s+5);
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);
format long
[POLES,ZEROS] = pzmap(sys)
format short
POLES = -5.000000000000007 -1.999999999999995 -1.000000000000001 -5.000000000000007 -1.999999999999995 -1.000000000000001 ZEROS = 0.000000000000000 + 1.414213562373095i 0.000000000000000 - 1.414213562373095i
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;
% _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 : 2.82e+04 ITE PFEAS DFEAS GFEAS PRSTATUS POBJ DOBJ MU TIME 0 1.0e+00 6.0e+00 6.2e+01 0.00e+00 -2.520000000e+02 0.000000000e+00 1.0e+00 0.00 1 1.6e-01 9.4e-01 2.4e+00 -1.22e+00 -5.799509574e+02 0.000000000e+00 1.6e-01 0.01 2 3.3e-02 1.9e-01 4.9e-01 -1.31e+00 -2.253278608e+02 0.000000000e+00 3.2e-02 0.01 3 4.5e-04 2.6e-03 1.2e-01 8.18e-01 -2.484152625e+00 0.000000000e+00 4.3e-04 0.01 4 2.0e-09 1.1e-08 1.3e-04 9.90e-01 -1.347953556e-05 0.000000000e+00 1.9e-09 0.01 5 1.9e-17 2.2e-15 1.1e-15 1.00e+00 -1.578114115e-13 0.000000000e+00 1.8e-17 0.01 Optimizer terminated. Time: 0.01 Interior-point solution summary Problem status : PRIMAL_AND_DUAL_FEASIBLE Solution status : OPTIMAL Primal. obj: -1.5781141147e-13 nrm: 7e-15 Viol. con: 3e-14 barvar: 0e+00 Dual. obj: 0.0000000000e+00 nrm: 2e+02 Viol. con: 0e+00 barvar: 1e-13 Optimizer summary Optimizer - time: 0.01 Interior-point - iterations : 5 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.2226 solvertime: 0.0198 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) = 806 + 0, nnz(ADA) = 9604, nnz(L) = 4851 it : b*y gap delta rate t/tP* t/tD* feas cg cg prec 0 : 6.08E+01 0.000 1 : 0.00E+00 1.76E+01 0.000 0.2899 0.9000 0.9000 2.00 1 1 1.4E+01 2 : 0.00E+00 4.37E+00 0.000 0.2479 0.9000 0.9000 1.25 1 1 3.6E+00 3 : 0.00E+00 1.02E+00 0.000 0.2336 0.9000 0.9000 1.11 1 1 1.2E+00 4 : 0.00E+00 3.67E-01 0.000 0.3588 0.9000 0.9000 1.01 1 1 8.2E-01 5 : 0.00E+00 1.39E-01 0.000 0.3801 0.9000 0.9000 0.97 1 1 6.8E-01 6 : 0.00E+00 4.80E-02 0.000 0.3447 0.9000 0.9000 0.98 1 1 6.3E-01 7 : 0.00E+00 1.67E-02 0.000 0.3479 0.9000 0.9000 0.99 1 1 4.6E-01 8 : 0.00E+00 6.14E-03 0.000 0.3676 0.9000 0.9000 0.99 1 1 1.7E-01 9 : 0.00E+00 2.16E-03 0.000 0.3524 0.9000 0.9000 0.99 1 1 6.2E-02 10 : 0.00E+00 8.23E-04 0.000 0.3801 0.9000 0.9000 0.99 1 1 2.4E-02 11 : 0.00E+00 2.77E-04 0.000 0.3362 0.9000 0.9000 0.99 1 1 8.2E-03 12 : 0.00E+00 9.95E-05 0.000 0.3596 0.9000 0.9000 0.98 1 1 3.0E-03 13 : 0.00E+00 3.21E-05 0.000 0.3232 0.9000 0.9000 0.96 1 1 1.0E-03 14 : 0.00E+00 1.14E-05 0.000 0.3537 0.9000 0.9000 0.94 1 1 3.8E-04 15 : 0.00E+00 3.83E-06 0.000 0.3369 0.9000 0.9000 0.89 1 1 1.4E-04 16 : 0.00E+00 1.42E-06 0.000 0.3714 0.9000 0.9000 0.82 1 1 6.2E-05 17 : 0.00E+00 4.92E-07 0.000 0.3457 0.9000 0.9000 0.70 1 1 2.8E-05 18 : 0.00E+00 1.90E-07 0.000 0.3864 0.9000 0.9000 0.53 1 1 1.6E-05 19 : 0.00E+00 6.47E-08 0.000 0.3407 0.9000 0.9000 0.38 1 1 7.1E-06 20 : 0.00E+00 2.48E-08 0.000 0.3829 0.9000 0.9000 0.22 1 1 1.1E-06 21 : 0.00E+00 8.40E-09 0.000 0.3389 0.9000 0.9000 0.15 1 1 3.8E-08 22 : 0.00E+00 3.18E-09 0.000 0.3785 0.9000 0.9000 0.08 1 1 2.4E-08 23 : 0.00E+00 1.08E-09 0.000 0.3397 0.9000 0.9000 0.06 1 1 1.3E-08 24 : 0.00E+00 4.13E-10 0.000 0.3821 0.9000 0.9000 0.02 1 1 8.3E-09 25 : 0.00E+00 1.39E-10 0.000 0.3367 0.9000 0.9000 0.02 1 1 4.6E-09 26 : 0.00E+00 5.23E-11 0.000 0.3760 0.9000 0.9000 0.00 1 1 2.9E-09 27 : 0.00E+00 1.76E-11 0.000 0.3370 0.9000 0.9000 0.01 1 1 1.6E-09 28 : 0.00E+00 6.45E-12 0.000 0.3660 0.9000 0.9000 -0.00 1 1 1.0E-09 29 : 0.00E+00 2.21E-12 0.000 0.3426 0.9000 0.9000 0.01 1 1 5.7E-10 iter seconds digits c*x b*y 29 0.3 11.7 -4.2054140828e-09 0.0000000000e+00 |Ax-b| = 5.1e-10, [Ay-c]_+ = 1.5E-11, |x|= 2.1e+03, |y|= 2.8e+03 Detailed timing (sec) Pre IPM Post 8.358E-03 1.264E-01 2.122E-03 Max-norms: ||b||=0, ||c|| = 1, Cholesky |add|=0, |skip| = 2, ||L.L|| = 2.54633. ans = struct with fields: yalmiptime: 0.1283 solvertime: 0.1372 info: 'Numerical problems (SeDuMi-1.3)' problem: 4 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ | ID| Constraint| Primal residual| Dual residual| ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ | #1| Matrix inequality| -3.9039e-12| 1.6778e-12| | #2| Matrix inequality| -1.5384e-11| 1.1196e-12| | #3| Matrix inequality| 5.1482e-05| 3.3443e-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.2234 0.1321 0.1050 0.0117 0.0405 0.0214 0.1321 0.4667 0.1365 0.0082 0.0293 0.0260 0.1050 0.1365 0.3740 0.0304 0.0425 0.0797 0.0117 0.0082 0.0304 0.2219 0.1318 0.1030 0.0405 0.0293 0.0425 0.1318 0.5661 0.1437 0.0214 0.0260 0.0797 0.1030 0.1437 0.3874 S = 231.6509 271.4625 138.7454 -12.7691 -26.6229 -28.3380 271.4625 695.2582 302.9054 -66.0838 49.7316 -82.3231 138.7454 302.9054 755.8996 -107.5396 -54.5168 148.4687 -12.7691 -66.0838 -107.5396 139.2652 37.6638 -0.3484 -26.6229 49.7316 -54.5168 37.6638 372.4472 4.0571 -28.3380 -82.3231 148.4687 -0.3484 4.0571 349.2328\begin{align} {\LARGE(1) \quad} L(\rho) = \left(\begin{array}{cc} 0.426\rho +3.28 & -4.62\rho -1.18 \\ -0.0386\rho -3.44 & 2.15\rho +1.71 \\ -0.201\rho -0.015 & 0.369\rho +2.13 \\ -0.528\rho -5.1 & 2.81\rho +13.6 \\ -0.0563\rho -1.28 & -0.931\rho -0.619 \\ 0.144\rho -1.43 & -0.0105\rho -0.579 \\ \end{array}\right) \end{align} \begin{align} {\LARGE(2) \quad} C_1(\rho) = \left(\begin{array}{cccccc} 5.5810^{-5}\rho +1.1210^{-4} & 1.7610^{-5}\rho +3.5110^{-5} & 1.4910^{-5}\rho +2.9810^{-5} & 2.9210^{-6}\rho +5.8310^{-6} & -1.6710^{-5}\rho -3.3410^{-5} & 1.7610^{-6}\rho +3.5210^{-6} \\ 2.9210^{-6}\rho +5.8310^{-6} & -3.0610^{-6}\rho -6.1210^{-6} & -1.3610^{-5}\rho -2.7310^{-5} & 5.5510^{-5}\rho +1.1110^{-4} & 6.5910^{-6}\rho +1.3210^{-5} & 7.0610^{-6}\rho +1.4110^{-5} \\ \end{array}\right) \end{align} \begin{align} {\LARGE(3) \quad} D_1(\rho) = \left(\begin{array}{cc} -2.2710^{-5}\rho -4.5410^{-5} & 5.3610^{-5}\rho +1.0710^{-4} \\ -4.2510^{-5}\rho -8.510^{-5} & 5.2710^{-5}\rho +1.0510^{-4} \\ \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 -2.7545 + 4.9849i -2.7545 - 4.9849i -3.9365 + 2.3963i -3.9365 - 2.3963i -0.9997 + 0.8124i -0.9997 - 0.8124i -1.4784 + 1.7608i -1.4784 - 1.7608i -0.8708 + 1.3187i -0.8708 - 1.3187i │ [ OK ] Zeros are all negative zeros if rho = -0.75 -3.0370 + 4.9222i -3.0370 - 4.9222i -3.6659 + 2.3590i -3.6659 - 2.3590i -1.0708 + 0.7836i -1.0708 - 0.7836i -1.4496 + 1.7451i -1.4496 - 1.7451i -0.9010 + 1.3292i -0.9010 - 1.3292i │ [ OK ] Zeros are all negative zeros if rho = -0.5 -3.3167 + 5.0392i -3.3167 - 5.0392i -3.3708 + 2.0407i -3.3708 - 2.0407i -1.1692 + 0.7549i -1.1692 - 0.7549i -1.4266 + 1.7320i -1.4266 - 1.7320i -0.9255 + 1.3353i -0.9255 - 1.3353i │ [ OK ] Zeros are all negative zeros if rho = -0.25 -3.5118 + 5.2765i -3.5118 - 5.2765i -3.0799 + 1.3499i -3.0799 - 1.3499i -1.3480 + 0.7328i -1.3480 - 0.7328i -1.4110 + 1.7214i -1.4110 - 1.7214i -0.9425 + 1.3369i -0.9425 - 1.3369i │ [ OK ] Zeros are all negative zeros if rho = 0 -3.6330 + 5.5611i -3.6330 - 5.5611i -4.0309 + 0.0000i -1.5197 + 0.0000i -1.6144 + 1.1227i -1.6144 - 1.1227i -1.4042 + 1.7134i -1.4042 - 1.7134i -0.9507 + 1.3341i -0.9507 - 1.3341i │ [ OK ] Zeros are all negative zeros if rho = 0.25 -3.7092 + 5.8651i -3.7092 - 5.8651i -4.8675 + 0.0000i -1.1911 + 0.0000i -1.3673 + 1.3387i -1.3673 - 1.3387i -1.4071 + 1.7078i -1.4071 - 1.7078i -0.9492 + 1.3270i -0.9492 - 1.3270i │ [ OK ] Zeros are all negative zeros if rho = 0.5 -3.7586 + 6.1791i -3.7586 - 6.1791i -5.4252 + 0.0000i -1.0614 + 0.0000i -1.1869 + 1.3845i -1.1869 - 1.3845i -1.4192 + 1.7044i -1.4192 - 1.7044i -0.9385 + 1.3159i -0.9385 - 1.3159i │ [ OK ] Zeros are all negative zeros if rho = 0.75 -3.7916 + 6.4994i -3.7916 - 6.4994i -5.8738 + 0.0000i -0.9851 + 0.0000i -1.4397 + 1.7031i -1.4397 - 1.7031i -1.0507 + 1.3787i -1.0507 - 1.3787i -0.9194 + 1.3009i -0.9194 - 1.3009i │ [ OK ] Zeros are all negative zeros if rho = 1 -3.8144 + 6.8240i -3.8144 - 6.8240i -6.2563 + 0.0000i -0.9361 + 0.0000i -1.4669 + 1.7039i -1.4669 - 1.7039i -0.9442 + 1.3499i -0.9442 - 1.3499i -0.8936 + 1.2821i -0.8936 - 1.2821i │ [ 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-4;
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(4) \quad}
\bar A(\rho) =
\left(\begin{array}{cc|cccccccccc}
-0.5\rho -5.6 & -0.085\rho -0.58 & 0.43\rho +1.910^{-3} & 0.1\rho +0.76 & -0.02\rho -0.38 & 4.410^{-3}\rho -0.062 & 0.11\rho -0.64 & 1.3\rho -1.8 & 0.18\rho +1.1 & -0.18\rho -1.1 & 0.97\rho -0.99 & 0.36\rho -1.2 \\
4.010^{-3}\rho -0.024 & 0.44\rho -5.7 & 9.710^{-3}\rho -0.15 & 0.029\rho -0.38 & -0.5\rho -0.011 & 0.83-0.14\rho & 0.19-0.031\rho & 0.56\rho +2.9 & 0.071\rho +2.5 & 0.19\rho +1.1 & 1.1\rho +5.9 & 0.11\rho -0.76 \\ \hline
6.8-1.1\rho & 0.017\rho -0.11 & 0.39\rho -2.4 & 0.3\rho -1.9 & 0.6-0.099\rho & 0.29-0.047\rho & 0.66\rho -4.0 & 0.21\rho -1.3 & 0.18\rho -1.1 & 1.910^{-16}\rho -1.210^{-15} & 1.2-0.2\rho & 0.1-0.016\rho \\
0 & 0 & 2.0-0.14\rho & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
3.810^{-17}\rho +2.310^{-16} & 1.1\rho +6.5 & 4.010^{-3}\rho +0.024 & 0.077\rho +0.45 & -0.4\rho -2.3 & -0.31\rho -1.8 & 5.110^{-16}\rho +3.010^{-15} & -0.049\rho -0.29 & -0.18\rho -1.0 & 0.68\rho +4.0 & 0.092\rho +0.54 & -0.086\rho -0.5 \\
0 & 0 & 0 & 0 & 0.058\rho +2.0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & -0.12\rho -8.0 & 2.6\rho -5.3 & 0.16\rho +0.39 & 0 & 2.3\rho +0.59 & 0.84\rho -2.2 \\
0 & 0 & 0 & 0 & 0 & 0 & 4.0-0.66\rho & 0.87-1.1\rho & -0.019\rho -1.7 & 0 & -1.1\rho -0.86 & 2.2-0.51\rho \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0.94-0.22\rho & -0.1\rho -7.510^{-3} & 0 & -0.18\rho -1.1 & 0.058\rho -0.52 \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & -1.1\rho -4.3 & -0.26\rho -2.5 & 0.038\rho -8.0 & -2.0\rho -11.0 & 1.7-0.28\rho \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0.49\rho +0.95 & -0.028\rho -0.64 & 0.68\rho +4.0 & 0.47\rho +0.31 & 0.27\rho +1.1 \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0.067\rho -1.0 & 0.71-0.072\rho & 0 & -0.064\rho -2.3 & 0.11\rho -1.2 \\
\end{array}\right)
\end{align}
\begin{align} {\LARGE(5) \quad}
\bar B(\rho) =
\left(\begin{array}{cc}
0.15\rho +0.3 & -2.410^{-3}\rho -4.810^{-3} \\
-5.110^{-18}\rho -1.010^{-17} & -0.15\rho -0.31 \\ \hline
0 & 0 \\
0 & 0 \\
0 & 0 \\
0 & 0 \\
0 & 0 \\
0 & 0 \\
0 & 0 \\
0 & 0 \\
0 & 0 \\
0 & 0 \\
\end{array}\right)
\end{align}
\begin{align} {\LARGE(6) \quad}
\bar C(\rho) =
\left(\begin{array}{cc|cccccccccc}
9.510^{-5}\rho +1.910^{-4} & -6.210^{-6}\rho -1.210^{-5} & -3.810^{-11}\rho -1.210^{-20} & 4.610^{-20}-2.110^{-12}\rho & 4.710^{-20}-3.110^{-12}\rho & -2.210^{-12}\rho -2.410^{-20} & -1.410^{-20}\rho -2.710^{-20} & -2.710^{-9}\rho -6.810^{-21} & 1.810^{-9}\rho +1.410^{-20} & -1.710^{-20}\rho -3.610^{-20} & 1.410^{-20}-7.710^{-10}\rho & -3.210^{-9}\rho -7.410^{-20} \\
4.910^{-6}\rho +9.910^{-6} & -9.110^{-5}\rho -1.810^{-4} & 3.010^{-12}\rho -8.810^{-20} & 1.410^{-11}\rho -7.210^{-20} & 2.710^{-20}-6.510^{-12}\rho & 1.410^{-20}-1.310^{-11}\rho & -4.210^{-20}\rho -8.610^{-20} & 5.010^{-20}-5.710^{-10}\rho & -8.510^{-10}\rho -4.110^{-20} & -2.710^{-20}\rho -4.110^{-20} & -1.410^{-9}\rho & 5.710^{-10}\rho -2.010^{-20} \\
\end{array}\right)
\end{align}
format long
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
format short
ans = 1.0e+03 * -2.210068849087277 -1.125268604692555 -0.705424243015968 -0.323071018438656 -0.465263208868048 -0.465185747360868 -0.000002089401602 -0.000001440227330 -0.000000291192170 -0.000000155926903 ans = 1.0e+03 * -2.324759230514694 -1.971755510337330 -0.357791331394135 -0.448895282616188 -0.465158906840635 -0.464427970407849 -0.000002531129411 -0.000001246362575 -0.000000153381884 -0.000000244542286
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 = linspace(-1,1,5)
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 = -1. │ [ OK ] Zero dynamics is stable if rho = -1, drho = -0.5. │ [ OK ] Zero dynamics is stable if rho = -1, drho = 0. │ [ OK ] Zero dynamics is stable if rho = -1, drho = 0.5. │ [ OK ] Zero dynamics is stable if rho = -1, drho = 1. │ [ OK ] Zero dynamics is stable if rho = -0.75, drho = -1. │ [ OK ] Zero dynamics is stable if rho = -0.75, drho = -0.5. │ [ OK ] Zero dynamics is stable if rho = -0.75, drho = 0. │ [ OK ] Zero dynamics is stable if rho = -0.75, drho = 0.5. │ [ OK ] Zero dynamics is stable if rho = -0.75, drho = 1. │ [ OK ] Zero dynamics is stable if rho = -0.5, drho = -1. │ [ OK ] Zero dynamics is stable if rho = -0.5, drho = -0.5. │ [ OK ] Zero dynamics is stable if rho = -0.5, drho = 0. │ [ OK ] Zero dynamics is stable if rho = -0.5, drho = 0.5. │ [ OK ] Zero dynamics is stable if rho = -0.5, drho = 1. │ [ OK ] Zero dynamics is stable if rho = -0.25, drho = -1. │ [ OK ] Zero dynamics is stable if rho = -0.25, drho = -0.5. │ [ OK ] Zero dynamics is stable if rho = -0.25, drho = 0. │ [ OK ] Zero dynamics is stable if rho = -0.25, drho = 0.5. │ [ OK ] Zero dynamics is stable if rho = -0.25, drho = 1. │ [ OK ] Zero dynamics is stable if rho = 0, drho = -1. │ [ OK ] Zero dynamics is stable if rho = 0, drho = -0.5. │ [ OK ] Zero dynamics is stable if rho = 0, drho = 0. │ [ OK ] Zero dynamics is stable if rho = 0, drho = 0.5. │ [ OK ] Zero dynamics is stable if rho = 0, drho = 1. │ [ OK ] Zero dynamics is stable if rho = 0.25, drho = -1. │ [ OK ] Zero dynamics is stable if rho = 0.25, drho = -0.5. │ [ OK ] Zero dynamics is stable if rho = 0.25, drho = 0. │ [ OK ] Zero dynamics is stable if rho = 0.25, drho = 0.5. │ [ OK ] Zero dynamics is stable if rho = 0.25, drho = 1. │ [ OK ] Zero dynamics is stable if rho = 0.5, drho = -1. │ [ OK ] Zero dynamics is stable if rho = 0.5, drho = -0.5. │ [ OK ] Zero dynamics is stable if rho = 0.5, drho = 0. │ [ OK ] Zero dynamics is stable if rho = 0.5, drho = 0.5. │ [ OK ] Zero dynamics is stable if rho = 0.5, drho = 1. │ [ OK ] Zero dynamics is stable if rho = 0.75, drho = -1. │ [ OK ] Zero dynamics is stable if rho = 0.75, drho = -0.5. │ [ OK ] Zero dynamics is stable if rho = 0.75, drho = 0. │ [ OK ] Zero dynamics is stable if rho = 0.75, drho = 0.5. │ [ OK ] Zero dynamics is stable if rho = 0.75, drho = 1. │ [ OK ] Zero dynamics is stable if rho = 1, drho = -1. │ [ OK ] Zero dynamics is stable if rho = 1, drho = -0.5. │ [ OK ] Zero dynamics is stable if rho = 1, drho = 0. │ [ OK ] Zero dynamics is stable if rho = 1, drho = 0.5. │ [ OK ] Zero dynamics is stable if rho = 1, drho = 1. │ [ 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 1 ans = logical 1