Tartalomjegyzék

LPV output passivation

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;

Instabil, nem min. fazisu MIMO LTI

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

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)/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
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                  : 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  8.8e-01  1.8e+00  -1.31e+00  -8.372130494e+02  0.000000000e+00   1.5e-01  0.01  
2   5.4e-02  1.4e-01  3.3e-01  -1.35e+00  -2.163114617e+02  0.000000000e+00   2.3e-02  0.01  
3   4.3e-04  1.1e-03  5.5e-02  8.66e-01   -1.369345252e+00  0.000000000e+00   1.8e-04  0.01  
4   3.5e-13  8.9e-13  9.2e-12  9.96e-01   -1.630210353e-09  0.000000000e+00   1.5e-13  0.01  
Optimizer terminated. Time: 0.02    


Interior-point solution summary
  Problem status  : PRIMAL_AND_DUAL_FEASIBLE
  Solution status : OPTIMAL
  Primal.  obj: -1.6302103533e-09   nrm: 4e-11    Viol.  con: 2e-10    barvar: 9e-18  
  Dual.    obj: 0.0000000000e+00    nrm: 1e+02    Viol.  con: 0e+00    barvar: 5e-11  
Optimizer summary
  Optimizer                 -                        time: 0.02    
    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.2242
    solvertime: 0.0250
          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) = 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.21E+00 0.000 0.2632 0.9000 0.9000   2.00  1  1  8.9E+00
  2 :   0.00E+00 2.82E+00 0.000 0.3430 0.9000 0.9000   1.26  1  1  3.2E+00
  3 :   0.00E+00 1.21E+00 0.000 0.4306 0.9000 0.9000   1.08  1  1  1.7E+00
  4 :   0.00E+00 5.93E-01 0.000 0.4893 0.9000 0.9000   1.04  1  1  1.1E+00
  5 :   0.00E+00 2.05E-01 0.000 0.3448 0.9000 0.9000   1.01  1  1  7.6E-01
  6 :   0.00E+00 7.53E-02 0.000 0.3681 0.9000 0.9000   1.00  1  1  6.5E-01
  7 :   0.00E+00 2.89E-02 0.000 0.3838 0.9000 0.9000   0.98  1  1  6.2E-01
  8 :   0.00E+00 9.92E-03 0.000 0.3432 0.9000 0.9000   0.97  1  1  5.3E-01
  9 :   0.00E+00 3.57E-03 0.000 0.3601 0.9000 0.9000   0.96  1  1  2.0E-01
 10 :   0.00E+00 1.21E-03 0.000 0.3379 0.9000 0.9000   0.95  1  1  7.2E-02
 11 :   0.00E+00 4.44E-04 0.000 0.3679 0.9000 0.9000   0.92  1  1  2.9E-02
 12 :   0.00E+00 1.49E-04 0.000 0.3346 0.9000 0.9000   0.88  1  1  1.1E-02
 13 :   0.00E+00 5.52E-05 0.000 0.3718 0.9000 0.9000   0.79  1  1  4.9E-03
 14 :   0.00E+00 1.89E-05 0.000 0.3427 0.9000 0.9000   0.67  1  1  2.2E-03
 15 :   0.00E+00 7.23E-06 0.000 0.3818 0.9000 0.9000   0.49  1  1  1.2E-03
 16 :   0.00E+00 2.49E-06 0.000 0.3452 0.9000 0.9000   0.35  1  1  5.3E-04
 17 :   0.00E+00 9.46E-07 0.000 0.3793 0.9000 0.9000   0.19  1  1  4.1E-06
 18 :   0.00E+00 3.22E-07 0.000 0.3407 0.9000 0.9000   0.14  1  1  2.2E-06
 19 :   0.00E+00 1.19E-07 0.000 0.3687 0.9000 0.9000   0.05  1  1  1.4E-06
 20 :   0.00E+00 4.04E-08 0.000 0.3397 0.9000 0.9000   0.06  1  1  7.4E-07
 21 :   0.00E+00 1.48E-08 0.000 0.3665 0.9000 0.9000   0.00  1  1  4.7E-07
 22 :   0.00E+00 5.06E-09 0.000 0.3420 0.9000 0.9000   0.04  1  1  2.6E-07
 23 :   0.00E+00 1.86E-09 0.000 0.3680 0.9000 0.9000  -0.02  1  1  1.6E-07
 24 :   0.00E+00 6.37E-10 0.000 0.3420 0.9000 0.9000   0.03  1  1  9.1E-08
 25 :   0.00E+00 2.34E-10 0.000 0.3674 0.9000 0.9000  -0.02  1  1  5.8E-08
 26 :   0.00E+00 7.98E-11 0.000 0.3409 0.9000 0.9000   0.02  1  1  3.2E-08
 27 :   0.00E+00 2.92E-11 0.000 0.3663 0.9000 0.9000  -0.03  1  2  2.1E-08
 28 :   0.00E+00 9.96E-12 0.000 0.3406 0.9000 0.9000   0.02  1  2  1.1E-08
 29 :   0.00E+00 3.65E-12 0.000 0.3662 0.9000 0.9000  -0.03  2  2  7.3E-09
 30 :   0.00E+00 1.24E-12 0.000 0.3409 0.9000 0.9000   0.02  2  2  4.0E-09
 31 :   0.00E+00 4.56E-13 0.000 0.3666 0.9000 0.9000  -0.03  2  2  2.6E-09
 32 :   0.00E+00 1.55E-13 0.000 0.3410 0.9000 0.9000   0.02  2  2  1.4E-09
 33 :   0.00E+00 5.70E-14 0.000 0.3665 0.9000 0.9000  -0.03  2  2  9.1E-10

iter seconds digits       c*x               b*y
 33      0.4  11.4 -4.3643988485e-07  0.0000000000e+00
|Ax-b| =   8.9e-10, [Ay-c]_+ =   3.1E-11, |x|=  1.1e+05, |y|=  3.1e+05

Detailed timing (sec)
   Pre          IPM          Post
8.067E-03    1.522E-01    2.444E-03    
Max-norms: ||b||=0, ||c|| = 1,
Cholesky |add|=0, |skip| = 2, ||L.L|| = 7.65576.
ans = 
  struct with fields:

    yalmiptime: 0.1251
    solvertime: 0.1635
          info: 'Numerical problems (SeDuMi-1.3)'
       problem: 4
 
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
|   ID|          Constraint|   Primal residual|   Dual residual|
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
|   #1|   Matrix inequality|       -4.0565e-11|      8.9323e-12|
|   #2|   Matrix inequality|       -3.4027e-11|      9.2439e-12|
|   #3|   Matrix inequality|        5.2031e-05|      1.5087e-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.2975    0.0494   -0.0952    0.0089    0.0946    0.0837
    0.0494    0.3665   -0.1271   -0.0580    0.1117    0.0901
   -0.0952   -0.1271    0.2795    0.0651   -0.0970   -0.1050
    0.0089   -0.0580    0.0651    0.5451    0.0279   -0.0237
    0.0946    0.1117   -0.0970    0.0279    0.3054    0.0232
    0.0837    0.0901   -0.1050   -0.0237    0.0232    0.3069
S =
   1.0e+05 *
    0.1098    0.0638   -0.0051    0.1726    0.0094   -0.1239
    0.0638    0.2504    0.0325    0.2648   -0.0724   -0.1064
   -0.0051    0.0325    0.1283    0.3389   -0.0283    0.0026
    0.1726    0.2648    0.3389    2.3854   -0.1701   -0.5792
    0.0094   -0.0724   -0.0283   -0.1701    0.0792   -0.0114
   -0.1239   -0.1064    0.0026   -0.5792   -0.0114    0.3070
\begin{align} {\LARGE(5) \quad} L(\rho) = \left(\begin{array}{cc} -4.72\rho -24.5 & 0.11\rho -7.53 \\ 0.957\rho -46.4 & 1.25\rho -10.9 \\ 2.36\rho +112.0 & 25.1-2.77\rho \\ -0.666\rho -45.9 & 1.2\rho -10.7 \\ 1.67\rho -126.0 & 3.33\rho -30.8 \\ -2.83\rho -116.0 & 2.89\rho -28.4 \\ \end{array}\right) \end{align} \begin{align} {\LARGE(6) \quad} C_1(\rho) = \left(\begin{array}{cccccc} 1.7410^{-4}-6.0410^{-5}\rho & 2.3510^{-5}\rho -1.8610^{-4} & 1.8710^{-4}\rho +3.5910^{-4} & -1.1210^{-6}\rho -1.9110^{-4} & 1.7310^{-4}\rho -3.7510^{-4} & 8.2810^{-6}\rho -1.6710^{-4} \\ -3.910^{-10}\rho -1.2110^{-4} & 1.5310^{-8}\rho +4.710^{-5} & 7.5510^{-9}\rho +3.7310^{-4} & 1.8710^{-8}\rho -2.2510^{-6} & 3.4610^{-4}-1.0310^{-8}\rho & 5.6110^{-10}\rho +1.6610^{-5} \\ \end{array}\right) \end{align} \begin{align} {\LARGE(7) \quad} D_1(\rho) = \left(\begin{array}{cc} 9.3310^{-5}\rho +1.5110^{-4} & 6.2810^{-5}\rho +1.1410^{-4} \\ 1.8710^{-4}-2.0210^{-8}\rho & 5.4110^{-9}\rho +1.2610^{-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
 -13.4718 + 0.0000i
  -4.0252 + 3.4690i
  -4.0252 - 3.4690i
  -3.4947 + 0.0000i
  -2.9024 + 0.0000i
  -0.8885 + 0.0000i
  -0.9692 + 1.1803i
  -0.9692 - 1.1803i
  -0.3061 + 0.4313i
  -0.3061 - 0.4313i
[  OK  ] Zeros are all negative

zeros if rho = -0.75
 -13.4125 + 0.0000i
  -3.9254 + 3.1196i
  -3.9254 - 3.1196i
  -3.3557 + 0.4947i
  -3.3557 - 0.4947i
  -0.8829 + 0.0000i
  -1.0098 + 1.1401i
  -1.0098 - 1.1401i
  -0.3013 + 0.4341i
  -0.3013 - 0.4341i
[  OK  ] Zeros are all negative

zeros if rho = -0.5
 -13.3005 + 0.0000i
  -3.6045 + 2.7449i
  -3.6045 - 2.7449i
  -3.7652 + 0.5596i
  -3.7652 - 0.5596i
  -0.8675 + 0.0000i
  -1.0502 + 1.0963i
  -1.0502 - 1.0963i
  -0.2967 + 0.4365i
  -0.2967 - 0.4365i
[  OK  ] Zeros are all negative

zeros if rho = -0.25
 -13.1260 + 0.0000i
  -5.1928 + 0.0000i
  -3.1104 + 2.6552i
  -3.1104 - 2.6552i
  -3.5702 + 0.0000i
  -0.8475 + 0.0000i
  -1.0904 + 1.0485i
  -1.0904 - 1.0485i
  -0.2922 + 0.4386i
  -0.2922 - 0.4386i
[  OK  ] Zeros are all negative

zeros if rho = 0
 -12.8703 + 0.0000i
  -6.4228 + 0.0000i
  -2.7690 + 2.7445i
  -2.7690 - 2.7445i
  -3.3492 + 0.0000i
  -0.8266 + 0.0000i
  -1.1305 + 0.9962i
  -1.1305 - 0.9962i
  -0.2879 + 0.4404i
  -0.2879 - 0.4404i
[  OK  ] Zeros are all negative

zeros if rho = 0.25
 -12.4941 + 0.0000i
  -7.4830 + 0.0000i
  -2.5334 + 2.8517i
  -2.5334 - 2.8517i
  -3.2056 + 0.0000i
  -0.8074 + 0.0000i
  -1.1703 + 0.9384i
  -1.1703 - 0.9384i
  -0.2838 + 0.4419i
  -0.2838 - 0.4419i
[  OK  ] Zeros are all negative

zeros if rho = 0.5
 -11.8821 + 0.0000i
  -8.6347 + 0.0000i
  -2.3559 + 2.9533i
  -2.3559 - 2.9533i
  -3.0869 + 0.0000i
  -0.7911 + 0.0000i
  -1.2100 + 0.8741i
  -1.2100 - 0.8741i
  -0.2799 + 0.4432i
  -0.2799 - 0.4432i
[  OK  ] Zeros are all negative

zeros if rho = 0.75
 -10.4831 + 0.8955i
 -10.4831 - 0.8955i
  -2.2157 + 3.0458i
  -2.2157 - 3.0458i
  -2.9804 + 0.0000i
  -0.7783 + 0.0000i
  -1.2494 + 0.8016i
  -1.2494 - 0.8016i
  -0.2762 + 0.4444i
  -0.2762 - 0.4444i
[  OK  ] Zeros are all negative

zeros if rho = 1
 -10.6759 + 2.0415i
 -10.6759 - 2.0415i
  -2.1022 + 3.1291i
  -2.1022 - 3.1291i
  -2.8808 + 0.0000i
  -0.7692 + 0.0000i
  -1.2887 + 0.7184i
  -1.2887 - 0.7184i
  -0.2727 + 0.4455i
  -0.2727 - 0.4455i
[  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-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) ]';

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(8) \quad} \bar A(\rho) = \left(\begin{array}{cc|cccccccccc} 0.061\rho -1.7 & 0.23\rho +0.016 & 2.7-5.910^{-3}\rho & 0.026\rho -1.5 & 9.610^{-3}\rho -0.013 & 0.032-0.049\rho & 2.5\rho -44.0 & -4.1\rho -18.0 & -1.4\rho -18.0 & 0.59\rho -133.0 & 3.2\rho +3.5 & 3.7\rho -64.0 \\ 0.28-0.034\rho & -0.14\rho -4.3 & 0.26-0.038\rho & -0.015\rho -0.15 & 0.16\rho +1.1 & 0.054\rho +1.7 & 1.8-0.59\rho & -0.76\rho -0.25 & 0.048-0.59\rho & 3.3-2.4\rho & 0.45\rho +1.1 & 2.9-0.82\rho \\ \hline 0.28\rho +5.8 & 0.26\rho -0.035 & -0.19\rho -2.4 & 0.091\rho +1.4 & 0.013\rho +0.067 & -0.061\rho -0.098 & -0.016\rho -1.1 & 0.043\rho +1.2 & -0.23\rho -2.1 & 0.059\rho +1.2 & 7.110^{-3}\rho +2.4 & -0.046\rho -1.0 \\ 1.7-0.45\rho & 2.810^{-3}\rho -0.01 & 0.39\rho -1.9 & 0.16-0.13\rho & 0.02-5.310^{-3}\rho & 7.710^{-3}\rho -0.029 & 0.088\rho -0.33 & 0.35-0.092\rho & 0.17\rho -0.62 & 0.35-0.093\rho & 0.72-0.19\rho & 0.082\rho -0.31 \\ 0.16-0.19\rho & -0.049\rho -2.6 & 0.069-0.081\rho & -0.056\rho -0.058 & 0.043\rho -0.55 & 2.210^{-3}\rho +0.69 & 0.029\rho -0.42 & 0.18-0.036\rho & 0.093\rho +1.2 & 0.023-0.039\rho & 1.2-0.059\rho & 0.033\rho -0.081 \\ -0.024\rho -6.310^{-3} & 0.036\rho +0.1 & -0.01\rho -2.710^{-3} & 2.210^{-3}-5.410^{-3}\rho & 8.610^{-3}\rho -0.38 & -9.110^{-3}\rho -0.027 & 0.01\rho +0.016 & -6.910^{-3}\rho -7.110^{-3} & -7.810^{-3}\rho -0.045 & -4.810^{-3}\rho -9.010^{-4} & -0.025\rho -0.046 & 5.110^{-3}\rho +3.110^{-3} \\ 0 & 0 & 0 & 0 & 0 & 0 & 0.41\rho +4.0 & 2.1\rho -6.310^{-3} & 1.3\rho +2.4 & 4.5\rho +24.0 & 1.3-1.3\rho & 1.1\rho +14.0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 17.0-1.4\rho & 1.2\rho +5.7 & 0.31\rho +6.8 & 45.0-1.0\rho & -0.93\rho -1.4 & 23.0-1.6\rho \\ 0 & 0 & 0 & 0 & 0 & 0 & 1.5\rho -37.0 & -4.5\rho -17.0 & -1.8\rho -19.0 & -2.0\rho -111.0 & 3.3\rho +3.3 & 2.5\rho -54.0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 15.0-0.73\rho & 1.8\rho +6.1 & 0.7\rho +7.2 & 0.53\rho +45.0 & -1.3\rho -1.3 & 23.0-1.1\rho \\ 0 & 0 & 0 & 0 & 0 & 0 & 43.0-2.6\rho & 3.5\rho +15.0 & 0.73\rho +15.0 & 122.0-1.9\rho & -2.8\rho -4.8 & 65.0-3.9\rho \\ 0 & 0 & 0 & 0 & 0 & 0 & 1.6\rho -40.0 & -4.8\rho -14.0 & -2.0\rho -17.0 & -2.4\rho -111.0 & 3.6\rho +0.85 & 2.3\rho -60.0 \\ \end{array}\right) \end{align} \begin{align} {\LARGE(9) \quad} \bar B(\rho) = \left(\begin{array}{cc} 6.010^{-3}\rho +1.8 & 0.012 \\ 0.97\rho +0.11 & 1.9 \\ \hline 5.610^{-16} & 0 \\ -3.110^{-16}\rho & -6.310^{-16} \\ 5.610^{-17}\rho & 1.110^{-16} \\ -1.110^{-17} & 0 \\ 0 & 0 \\ 0 & 0 \\ 0 & 0 \\ 0 & 0 \\ 0 & 0 \\ 0 & 0 \\ \end{array}\right) \end{align} \begin{align} {\LARGE(10) \quad} \bar C(\rho) = \left(\begin{array}{cc|cccccccccc} 9.310^{-4}-2.110^{-6}\rho & 4.010^{-4}\rho +3.710^{-5} & 1.710^{-12}-2.110^{-11}\rho & 8.310^{-13}-1.010^{-11}\rho & 2.710^{-12}\rho -2.210^{-13} & 1.610^{-13}-2.010^{-12}\rho & 1.310^{-11}-1.610^{-10}\rho & 1.110^{-12}-1.310^{-11}\rho & 4.510^{-12}-5.710^{-11}\rho & 3.410^{-11}-4.310^{-10}\rho & 1.310^{-12}-1.610^{-11}\rho & 2.210^{-11}-2.710^{-10}\rho \\ -2.110^{-11}\rho -4.110^{-6} & 8.010^{-4}-5.810^{-12}\rho & 4.210^{-11}-6.310^{-12}\rho & 3.310^{-12}\rho +2.110^{-11} & 2.410^{-12}\rho -5.510^{-12} & 2.010^{-13}\rho +4.010^{-12} & 3.210^{-10}-3.910^{-10}\rho & 1.510^{-8}\rho +2.610^{-11} & 7.610^{-9}\rho +1.110^{-10} & 1.910^{-8}\rho +8.610^{-10} & 3.210^{-11}-1.010^{-8}\rho & 5.410^{-10}-5.610^{-10}\rho \\ \end{array}\right) \end{align}

Check stability of the transformed system

(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
Output:
ans =
   1.0e+04 *
   -3.6239
   -3.3649
   -0.7888
   -1.3679
   -2.1188
   -2.1188
   -0.0000
   -0.0000
   -0.0000
   -0.0000
ans =
   1.0e+04 *
   -3.4864
   -3.9769
   -0.9946
   -2.1395
   -2.1188
   -2.1190
   -0.0000
   -0.0000
   -0.0000
   -0.0000

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