Tartalomjegyzék

LPV output passivation

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;

Instabil, nem min. fazisu MIMO LTI

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
Output:
POLES =
  -5.000000000000007
  -1.999999999999995
  -1.000000000000001
  -5.000000000000007
  -1.999999999999995
  -1.000000000000001
ZEROS =
  0.000000000000000 + 1.414213562373095i
  0.000000000000000 - 1.414213562373095i

Csinalok belole LPV-t

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
Output:
n_x =
     6
n_u =
     2
n_y =
     2
n_r =
     2
n_yp =
     2

Visszacsatolasi erosites

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
Output:
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

Observer design

Változók és mátrixok deklarálása

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
    ]
Output:
++++++++++++++++++++++++++++++++++
|   ID|                Constraint|
++++++++++++++++++++++++++++++++++
|   #1|   Matrix inequality 16x16|
|   #2|   Matrix inequality 16x16|
|   #3|   Matrix inequality 12x12|
++++++++++++++++++++++++++++++++++

Solve the optimization problem

sdpopts = sdpsettings('solver','sedumi');
optimize(Constraints,[],sdpopts)
check(Constraints)
Output:
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|
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 

Retrive values

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};
Output:
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}

Open loop matrices

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
Output:
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

Zero dynamics of the obtained LPV

Generate transformation matrix

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) ]';

Transformed LPV system

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}

Check stability of the transformed system

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
Output:
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

Symbolic zero dynamics

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')
Output:
│   [  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)) )
Output:
is_negdef =
  function_handle with value:
    @(A)all(round(eig(A),prec)<=0)
ans =
  logical
   1
ans =
  logical
   1