Teljes Matlab script kiegészítő függvényekkel.
File: 2018_01_25_K_prelim_L_codesign_v5_rhofuggoBL.m Directory: projects/3_outsel/2017_11_13_lpv_passivity Author: Peter Polcz (ppolcz@gmail.com)
Created on 2018. January 25.
Inherited from
LPV output passivization
File: d2018_01_09_K_prelim_L_codesign_v4.m Directory: projects/3_outsel/2017_11_13_lpv_passivity Author: Peter Polcz (ppolcz@gmail.com)
Created on 2018. January 09. Modified on 2018. January 17. Modified on 2018. January 25.
FIGYELEM, a rendszermodell, nem stabil!
s = tf('s');
% unstable MIMO
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\begin{align} {\LARGE(1) \quad} H(s) = \left(\begin{array}{cc} -\frac{s-1}{-s^2+s+2} & \frac{10}{10s^2+29s-3} \\ \frac{s-7}{s^2+6s+5} & \frac{s-6}{s^2+5s+6} \\ \end{array}\right) \end{align}
Az előző LTI modellt egy kicsit megperturbálom (kezdetben csak az A mátrixot).
$$ \begin{aligned} &\Sigma: \left\{\begin{aligned} &\dot x = A(\rho) x + B(\rho) u,~~~ \rho \in \mathcal P \\ &y = C x \\ \end{aligned}\right. \\ &\begin{aligned} \text{where: } & A(\rho) = A_0 + A_1 \rho \in \mathbb{R}^{n\times n} \\ & B(\rho) = B_0 + B_1 \rho \in \mathbb{R}^{n\times r} \\ & C \in \mathbb{R}^{m\times n} \\ & D = 0_{m\times r} \end{aligned} \end{aligned} $$
rho_lim = [
-1 1
];
A1 = A0;
A1(abs(A0) < 1) = 0;
A1 = A1 .* randn(size(A1))/10;
B1 = [ 2\B0(:,1)+3\B0(:,2) B0(:,1)*0 ]*0;
A_fh = @(rho) A0 + rho*A1;
B_fh = @(rho) B0 + rho*B1;
\begin{align} {\LARGE(2) \quad}
A_0 = \left(\begin{array}{cccccc}
-5.37 & -0.933 & -0.464 & 0.0928 & 0.464 & 0.464 \\
1.87 & -0.319 & -0.188 & 0.0375 & 0.188 & 0.188 \\
-0.292 & -0.731 & -1.35 & 0.29 & -1.55 & -1.55 \\
0.0585 & 0.146 & 0.224 & 0.0552 & 0.276 & 0.276 \\
0.292 & 0.731 & -3.19 & 0.638 & -1.81 & 0.188 \\
0.292 & 0.731 & 0.101 & -0.0203 & 1.9 & -0.101 \\
\end{array}\right)
,\quad
A_1 = \left(\begin{array}{cccccc}
-0.323 & 0 & 0 & 0 & 0 & 0 \\
0.111 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0.0044 & 0 & 0.154 & 0.236 \\
0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0.136 & 0 & 0.116 & 0 \\
0 & 0 & 0 & 0 & 0.343 & 0 \\
\end{array}\right)
\end{align}
\begin{align} {\LARGE(3) \quad}
B_0 = \left(\begin{array}{cc}
2 & 0 \\
0 & 0 \\
1.15 & 2 \\
-0.229 & 0 \\
-1.15 & 2 \\
-1.15 & 0 \\
\end{array}\right)
,\quad
B_1 = \left(\begin{array}{cc}
0 & 0 \\
0 & 0 \\
0 & 0 \\
0 & 0 \\
0 & 0 \\
0 & 0 \\
\end{array}\right)
\end{align}
\begin{align} {\LARGE(4) \quad}
C = \left(\begin{array}{cccccc}
0.169 & 0.422 & 0.256 & 0.949 & -0.256 & -0.256 \\
0.701 & -1.25 & -0.444 & 0.0889 & 0.944 & -1.06 \\
\end{array}\right)
,\quad
D = \left(\begin{array}{cc}
0 & 0 \\
0 & 0 \\
\end{array}\right)
\end{align}
Matrix $K$ is preliminarily design for the nominal value of $\rho$, than checked whether the closed loop is stable for other $\rho \in \mathcal P$ values. Let $G = I_m$.
G = eye(n_u);
K = place(A0,B0,linspace(-5,-6,n_x));
Free matrix variables. All matrices are assumed to be parameter independent.
C1 = sdpvar(n_yp,n_x,'full');
D1 = sdpvar(n_yp,n_y,'full');
Q = sdpvar(n_x,n_x,'symmetric');
S = sdpvar(n_x,n_x,'symmetric');
N_1 = sdpvar(n_x,n_y,'full');
N_2 = sdpvar(n_x,n_y,'full');
N_fh = @(rho) N_1 + rho*N_2;
P = blkdiag(Q,S);
where $N = S L$.
$$ \begin{align} &\wt A(\rho) = \pmqty{ A(\rho) & 0 \\ 0 & A(\rho) - L C },~ \wt B = \pmqty{B \\ 0} \\ &\wt C = \pmqty{ D_1 C + C_1 & -C_1 }. \end{align} $$
Ao = @(L,rho) [
A_fh(rho) zeros(n_x,n_x)
zeros(n_x,n_x) A_fh(rho)-L(rho)*C
];
Bo = @(rho) [B_fh(rho) ; zeros(n_x,n_u)];
Co = [D1*C+C1 -C1];
Do = zeros(n_yp,n_r);
$$ \begin{aligned} \wt A_c(\rho) = \begin{pmatrix} A(\rho)-B(\rho) K & B(\rho) K \\ 0 & A(\rho) - L C \end{pmatrix},~ \wt B_c(\rho) = \begin{pmatrix} B(\rho) G \\ 0 \end{pmatrix}. \end{aligned} $$
Ac = @(L,rho) [
A_fh(rho)-B_fh(rho)*K B_fh(rho)*K
zeros(n_x,n_x) A_fh(rho)-L(rho)*C
];
Bc = @(rho) Bo(rho)*G;
W = eye(n_yp);
$$ \begin{aligned} &\wt A_c({\color{red} \rho})^T P + P \wt A_c({\color{red} \rho}) = \spmqty{ {\blue Q} A_K({\red \rho}) + A_K^T({\red \rho}) {\blue Q} & {\blue Q} B({\red \rho}) K \\ K^T B^T({\red \rho}) {\blue Q} & {\blue S} A({\red \rho}) + A({\red \rho})^T {\blue S} - {\blue N}({\red \rho}) C - C^T {\blue N^T}({\red \rho}), } \\ &\text{where } A_K({\red \rho}) = A({\red \rho}) - B({\red \rho}) K. \nonumber \end{aligned} $$
AcP_PAc = @(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)'
];
$$ M_2 = \spmqty{ { \wt A_c({\color{red} \rho})^T} { P} + { P} { \wt A_c({\color{red} \rho})} & { P} \wt B_c({\red \rho}) - { \wt C^T({\red \rho})} & { \wt C^T({\red \rho})} \\ \wt B_c^T({\red \rho}) { P}- { \wt C({\red \rho})} & 0 & 0 \\ { \wt C({\red \rho})} & 0 & -W^{-1} } $$
M2 = @(rho) [
AcP_PAc(rho) P*Bc(rho)-Co' Co'
Bc(rho)'*P-Co zeros(n_r,n_r) zeros(n_r,n_yp)
Co zeros(n_yp,n_r) -inv(W)
];
% Constraints
Constraints = [
M2(rho_lim(1)) <= 0
M2(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| ++++++++++++++++++++++++++++++++++
optimize(Constraints)
check(Constraints)
Q = value(Q);
S = value(S);
N_1 = value(N_1);
N_2 = value(N_2);
L_1 = S\N_1;
L_2 = S\N_2;
L_fh = @(rho) L_1 + rho*L_2;
P = value(P);
Co = value(Co);
C1 = value(C1);
D1 = value(D1);
Problem Name : Objective sense : min Type : CONIC (conic optimization problem) Constraints : 82 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 : 82 Cones : 0 Scalar variables : 0 Matrix variables : 3 Integer variables : 0 Optimizer - threads : 4 Optimizer - solved problem : the primal Optimizer - Constraints : 82 Optimizer - Cones : 0 Optimizer - Scalar variables : 0 conic : 0 Optimizer - Semi-definite variables: 3 scalarized : 350 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 : 3403 after factor : 3403 Factor - dense dim. : 0 flops : 5.86e+05 ITE PFEAS DFEAS GFEAS PRSTATUS POBJ DOBJ MU TIME 0 2.4e+00 1.0e+00 5.0e+00 0.00e+00 3.998800000e+00 0.000000000e+00 1.0e+00 0.01 1 5.2e-01 2.2e-01 2.3e+00 1.00e+00 8.755693332e-01 0.000000000e+00 2.2e-01 0.02 2 1.9e-01 7.9e-02 1.4e+00 9.99e-01 3.166073696e-01 0.000000000e+00 7.9e-02 0.02 3 7.3e-02 3.1e-02 8.7e-01 9.96e-01 1.231920274e-01 0.000000000e+00 3.1e-02 0.02 4 2.7e-02 1.1e-02 5.3e-01 9.96e-01 4.570574136e-02 0.000000000e+00 1.1e-02 0.02 5 5.8e-03 2.5e-03 2.4e-01 9.92e-01 9.875092415e-03 0.000000000e+00 2.5e-03 0.02 6 1.2e-03 5.0e-04 1.1e-01 9.82e-01 2.040980985e-03 0.000000000e+00 5.0e-04 0.02 7 3.8e-04 1.6e-04 6.0e-02 9.67e-01 6.547027256e-04 0.000000000e+00 1.6e-04 0.02 8 1.1e-04 4.7e-05 3.2e-02 9.74e-01 1.972824594e-04 0.000000000e+00 4.7e-05 0.02 9 4.4e-05 1.9e-05 2.0e-02 9.66e-01 7.865094703e-05 0.000000000e+00 1.9e-05 0.02 10 9.3e-06 3.9e-06 8.7e-03 9.46e-01 1.700041537e-05 0.000000000e+00 3.9e-06 0.03 11 2.6e-06 1.1e-06 4.2e-03 8.88e-01 5.056939028e-06 0.000000000e+00 1.1e-06 0.03 12 7.5e-07 3.2e-07 1.8e-03 7.49e-01 1.585240685e-06 0.000000000e+00 3.2e-07 0.03 13 3.1e-07 1.3e-07 8.7e-04 5.06e-01 7.471267729e-07 0.000000000e+00 1.3e-07 0.03 14 9.2e-08 3.9e-08 3.6e-04 4.43e-01 2.384750889e-07 0.000000000e+00 3.9e-08 0.03 15 3.0e-08 1.3e-08 1.1e-04 1.94e-01 -1.371440162e-08 0.000000000e+00 1.3e-08 0.03 16 6.0e-09 2.5e-09 2.9e-05 2.22e-01 -8.385578824e-08 0.000000000e+00 2.5e-09 0.03 17 1.6e-09 6.6e-10 7.1e-06 6.03e-02 -1.577668375e-07 0.000000000e+00 6.6e-10 0.03 18 4.3e-10 1.8e-10 2.0e-06 1.23e-02 -1.889399310e-07 0.000000000e+00 1.8e-10 0.03 19 1.7e-10 7.2e-11 7.2e-07 -4.84e-02 -2.287530519e-07 0.000000000e+00 7.2e-11 0.04 20 4.7e-11 2.0e-11 2.4e-07 1.02e-01 -1.594187513e-07 0.000000000e+00 2.0e-11 0.04 21 1.3e-11 5.4e-12 5.3e-08 -4.32e-02 -2.478096800e-07 0.000000000e+00 5.4e-12 0.04 22 2.8e-12 1.2e-12 1.4e-08 6.76e-02 -1.809508529e-07 0.000000000e+00 1.2e-12 0.04 23 7.8e-13 3.3e-13 3.8e-09 3.36e-02 -1.877296559e-07 0.000000000e+00 3.3e-13 0.04 24 2.3e-13 9.8e-14 1.1e-09 1.21e-02 -1.940468187e-07 0.000000000e+00 9.7e-14 0.04 25 9.2e-14 4.2e-14 3.9e-10 -9.85e-02 -2.503239479e-07 0.000000000e+00 3.9e-14 0.04 26 2.4e-14 1.9e-14 1.3e-10 1.02e-01 -1.599584278e-07 0.000000000e+00 1.0e-14 0.04 27 8.3e-15 8.2e-15 3.7e-11 -1.45e-02 -2.216619486e-07 0.000000000e+00 3.4e-15 0.04 28 9.0e-15 3.6e-15 7.7e-12 -4.67e-02 -2.078673001e-07 0.000000000e+00 7.2e-16 0.04 29 9.6e-15 1.4e-14 2.4e-12 6.01e-02 -1.881283600e-07 0.000000000e+00 2.1e-16 0.05 30 3.5e-14 8.0e-15 7.0e-13 1.78e-02 -1.913183009e-07 0.000000000e+00 6.1e-17 0.05 31 1.3e-14 3.9e-15 2.7e-13 -2.11e-02 -2.022821419e-07 0.000000000e+00 2.4e-17 0.05 32 7.8e-15 4.2e-15 6.4e-14 -4.37e-02 -1.962618795e-07 0.000000000e+00 5.7e-18 0.05 33 3.1e-14 3.9e-15 2.8e-14 5.85e-02 -2.018576626e-07 0.000000000e+00 2.5e-18 0.05 34 8.6e-15 5.5e-15 6.0e-15 -5.72e-02 -2.061561207e-07 0.000000000e+00 5.7e-19 0.05 35 1.8e-14 5.2e-15 2.6e-15 4.81e-02 -2.135388074e-07 0.000000000e+00 2.5e-19 0.05 36 4.1e-15 3.4e-15 5.8e-16 -6.22e-02 -2.161701279e-07 0.000000000e+00 5.8e-20 0.05 37 1.7e-15 4.0e-15 2.4e-16 5.55e-02 -2.215534097e-07 0.000000000e+00 2.4e-20 0.06 38 3.8e-16 6.2e-15 5.4e-17 -5.73e-02 -2.221187378e-07 0.000000000e+00 5.7e-21 0.06 39 1.6e-16 3.5e-15 2.2e-17 5.83e-02 -2.345455988e-07 0.000000000e+00 2.4e-21 0.06 40 3.6e-17 8.9e-15 4.9e-18 -7.84e-02 -2.349898597e-07 0.000000000e+00 5.7e-22 0.06 41 1.6e-17 7.2e-15 2.2e-18 7.95e-02 -2.445369317e-07 0.000000000e+00 2.4e-22 0.06 42 3.5e-18 1.4e-14 4.7e-19 -7.27e-02 -2.471943501e-07 0.000000000e+00 5.8e-23 0.07 43 1.6e-18 4.9e-15 2.1e-19 7.82e-02 -2.580022299e-07 0.000000000e+00 2.5e-23 0.07 44 3.5e-19 6.7e-15 4.5e-20 -6.49e-02 -2.621010202e-07 0.000000000e+00 5.9e-24 0.07 45 1.5e-19 4.9e-15 1.9e-20 6.77e-02 -2.708044697e-07 0.000000000e+00 2.5e-24 0.07 46 3.3e-20 5.8e-15 4.2e-21 -5.55e-02 -2.739479729e-07 0.000000000e+00 6.0e-25 0.07 47 2.5e-17 7.1e-15 3.0e-21 5.69e-02 -2.779917531e-07 0.000000000e+00 4.3e-25 0.08 48 2.5e-17 5.8e-15 3.0e-21 1.59e-02 -2.780344052e-07 0.000000000e+00 4.2e-25 0.08 49 2.5e-17 8.4e-15 3.0e-21 1.54e-02 -2.780370327e-07 0.000000000e+00 4.2e-25 0.08 50 2.5e-17 1.1e-14 3.0e-21 1.55e-02 -2.780396983e-07 0.000000000e+00 4.2e-25 0.08 51 2.5e-17 1.1e-14 3.0e-21 1.61e-02 -2.780396983e-07 0.000000000e+00 4.2e-25 0.09 Optimizer terminated. Time: 0.10 Interior-point solution summary Problem status : PRIMAL_AND_DUAL_FEASIBLE Solution status : OPTIMAL Primal. obj: -2.7803969833e-07 nrm: 5e+09 Viol. con: 3e-07 barvar: 2e-07 Dual. obj: 0.0000000000e+00 nrm: 6e+09 Viol. con: 0e+00 barvar: 1e-05 Optimizer summary Optimizer - time: 0.10 Interior-point - iterations : 52 time: 0.09 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.3025 solvertime: 0.1106 info: 'Successfully solved (MOSEK)' problem: 0 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ | ID| Constraint| Primal residual| Dual residual| ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ | #1| Matrix inequality| -1.5857e-07| 1.424e-16| | #2| Matrix inequality| -1.1524e-07| 1.0903e-16| | #3| Matrix inequality| 2.4627e-05| 3.7307e-17| ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
[ OK ] M2(0) <= 0 [ OK ] M2(0.1) <= 0 [ OK ] M2(-0.1) <= 0 [ OK ] P > 0 [ OK ] A - LC < 0 [ OK ] norm(D1) > 1e-3 [ OK ] Zeros of the open loop system are stable [ OK ] Zeros of the closed loop system are stable
Written on 2018. January 17.
Ebben a részben $A(\rho)$, $B(\rho)$, $C$, $D$ új értelmet nyer. Ezek adják meg az observerrel kiegészített dinamika mátrixait.
A_fh = @(rho) Ao(L_fh,rho);
B_fh = Bo;
A0 = A_fh(0);
A1 = A_fh(1) - A0;
B0 = B_fh(0);
B1 = B_fh(1) - B0;
C0 = Co;
C1 = C0*0;
D0 = Do;
D1 = Do*0;
A_fh = @(rho) A0 + A1*rho;
B_fh = @(rho) B0 + B1*rho;
C_fh = @(rho) C0 + C1*rho;
D_fh = @(rho) D0 + D1*rho;
A = {A0, A1};
B = {B0, B1};
C = {C0, C1};
D = {D0, D1};
AA = [A{:}];
BB = [B{:}];
\begin{align} {\LARGE(5) \quad}
A(\rho) =
\left(\begin{array}{cccccccccccc}
-0.32\rho -5.4 & -0.93 & -0.46 & 0.093 & 0.46 & 0.46 & 0 & 0 & 0 & 0 & 0 & 0 \\
0.11\rho +1.9 & -0.32 & -0.19 & 0.038 & 0.19 & 0.19 & 0 & 0 & 0 & 0 & 0 & 0 \\
-0.29 & -0.73 & 4.410^{-3}\rho -1.4 & 0.29 & 0.15\rho -1.5 & 0.24\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.14\rho -3.2 & 0.64 & 0.12\rho -1.8 & 0.19 & 0 & 0 & 0 & 0 & 0 & 0 \\
0.29 & 0.73 & 0.1 & -0.02 & 0.34\rho +1.9 & -0.1 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 0.61\rho +2.2 & -0.89\rho -0.19 & 1.9-0.21\rho & 1.1\rho +19.0 & 0.74\rho +1.2 & -1.4\rho -11.0 \\
0 & 0 & 0 & 0 & 0 & 0 & 2.1\rho +9.9 & 6.7-1.2\rho & 5.4-0.11\rho & 3.2\rho +28.0 & 1.1\rho -3.2 & -3.0\rho -12.0 \\
0 & 0 & 0 & 0 & 0 & 0 & -5.1\rho -21.0 & 3.9\rho -20.0 & 0.64\rho -16.0 & -7.2\rho -73.0 & 7.9-3.2\rho & 7.9\rho +29.0 \\
0 & 0 & 0 & 0 & 0 & 0 & 2.1\rho +8.9 & 6.6-1.6\rho & 5.8-0.26\rho & 2.9\rho +30.0 & 1.4\rho -2.7 & -3.1\rho -13.0 \\
0 & 0 & 0 & 0 & 0 & 0 & 5.6\rho +25.0 & 16.0-3.3\rho & 11.0-0.089\rho & 9.2\rho +81.0 & 3.2\rho -8.2 & -8.4\rho -38.0 \\
0 & 0 & 0 & 0 & 0 & 0 & 5.1\rho +24.0 & 15.0-3.9\rho & 13.0-0.66\rho & 7.1\rho +75.0 & 3.8\rho -3.6 & -7.6\rho -36.0 \\
\end{array}\right)
\end{align}
\begin{align} {\LARGE(6) \quad}
B(\rho) =
\left(\begin{array}{cc}
2.0 & 0 \\
0 & 0 \\
1.1 & 2.0 \\
-0.23 & 0 \\
-1.1 & 2.0 \\
-1.1 & 0 \\
0 & 0 \\
0 & 0 \\
0 & 0 \\
0 & 0 \\
0 & 0 \\
0 & 0 \\
\end{array}\right)
\end{align}
\begin{align} {\LARGE(7) \quad}
C(\rho) = C = 10^{-3}
\left(\begin{array}{cccccccccccc}
0.19 & -0.16 & 0.36 & 1.3 & -0.19 & -0.49 & 0.2 & 0.33 & -0.18 & -0.11 & 0.13 & -0.097 \\
0.091 & 0.078 & 0.62 & 0.7 & 0.14 & 0.18 & 2.610^{-3} & 0.19 & -0.46 & -0.13 & -0.3 & -0.33 \\
\end{array}\right)
\end{align}
Nálam az $E_c$ most az $Im(B)$.
Im_B = IMA([B{1} B{2}]);
\begin{align} {\LARGE(8) \quad}
\mathrm{Im}\big(B(\rho)\big) = \left(\begin{array}{cc}
-0.707 & 0 \\
0 & 0 \\
-0.406 & -0.707 \\
0.0811 & 0 \\
0.406 & -0.707 \\
0.406 & 0 \\
0 & 0 \\
0 & 0 \\
0 & 0 \\
0 & 0 \\
0 & 0 \\
0 & 0 \\
\end{array}\right)
\end{align}
Ker_C = INTS(KER(C{1}), KER(C{2}));
\begin{align} {\LARGE(9) \quad}
\mathrm{Ker}(C) = \left(\begin{array}{cccccccccc}
-0.992 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
-0.0143 & 0.978 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0.028 & -0.0591 & -0.843 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0.105 & 0.0804 & 0.277 & 0.448 & 0 & 0 & 0 & 0 & 0 & 0 \\
-0.0165 & -0.06 & 0.115 & 0.13 & 0.941 & 0 & 0 & 0 & 0 & 0 \\
-0.0426 & -0.121 & 0.184 & 0.444 & -0.264 & 0.595 & 0 & 0 & 0 & 0 \\
0.017 & 0.0341 & -0.0252 & -0.223 & 0.0832 & 0.376 & -0.878 & 0 & 0 & 0 \\
0.0274 & 0.0201 & 0.0756 & -0.469 & 0.0794 & 0.467 & 0.356 & 0.611 & 0 & 0 \\
-0.0129 & 0.0599 & -0.268 & 0.442 & 0.0711 & 0.0554 & -0.04 & 0.401 & -0.616 & 0 \\
-0.0091 & 0.0069 & -0.0704 & 0.198 & -0.005 & -0.1 & -0.0926 & 0.323 & 0.451 & -0.786 \\
0.0129 & 0.0834 & -0.212 & 0.0148 & 0.153 & 0.514 & 0.304 & -0.571 & 0.0364 & -0.256 \\
-0.0068 & 0.0477 & -0.195 & 0.283 & 0.0628 & 0.094 & 0.0093 & 0.19 & 0.645 & 0.563 \\
\end{array}\right)
\end{align}
We compute $\mathcal V^*$ and $\mathcal R^*$.
[R,V] = CSA([A{1} A{2}], Im_B, Ker_C);
dim_V = size(V,2);
dim_Ker_V = size(V,1) - size(V,2);
assert(rank(R) == 0 && rank([V Ker_C]) == rank(Ker_C), ...
'Strong invertibility condition is not satisfied!');
assert(rank(INTS(V, Im_B)) == 0)
We ensure that $\mathcal R^* = \{0\}$ and that $\mathcal V^* \subseteq \mathrm{Ker}(C)$. Furthermore, we checked the strong invertibility condition $\mathcal V^* \cap \mathrm{Im}\big(B(\rho)\big) = \{0\}$.
\begin{align} {\LARGE(10) \quad} \mathcal V^* = \left(\begin{array}{cccccccccc} -0.992 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ -0.0143 & 0.978 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0.028 & -0.0591 & -0.843 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0.105 & 0.0804 & 0.277 & 0.448 & 0 & 0 & 0 & 0 & 0 & 0 \\ -0.0165 & -0.06 & 0.115 & 0.13 & 0.941 & 0 & 0 & 0 & 0 & 0 \\ -0.0426 & -0.121 & 0.184 & 0.444 & -0.264 & 0.595 & 0 & 0 & 0 & 0 \\ 0.017 & 0.0341 & -0.0252 & -0.223 & 0.0832 & 0.376 & -0.878 & 0 & 0 & 0 \\ 0.0274 & 0.0201 & 0.0756 & -0.469 & 0.0794 & 0.467 & 0.356 & 0.611 & 0 & 0 \\ -0.0129 & 0.0599 & -0.268 & 0.442 & 0.0711 & 0.0554 & -0.04 & 0.401 & -0.616 & 0 \\ -0.0091 & 0.0069 & -0.0704 & 0.198 & -0.005 & -0.1 & -0.0926 & 0.323 & 0.451 & -0.786 \\ 0.0129 & 0.0834 & -0.212 & 0.0148 & 0.153 & 0.514 & 0.304 & -0.571 & 0.0364 & -0.256 \\ -0.0068 & 0.0477 & -0.195 & 0.283 & 0.0628 & 0.094 & 0.0093 & 0.19 & 0.645 & 0.563 \\ \end{array}\right) \end{align}1. variáns
Let $T = \begin{pmatrix} {V^*}^\perp & L \end{pmatrix}^T$, where ${V^*}^\perp$ and $L$ are orthonormal bases for ${\mathcal V^*}^\perp \not\perp \mathcal L \subseteq \mathrm{Im}\big(B(\rho)\big)$, respectively.
L = ORTCO(Im_B);
dim_L = size(L,2);
assert(dim_L >= dim_V);
T = [ ORTCO(V) L(:,1:dim_V) ]';
\begin{align} {\LARGE(11) \quad}
{\mathcal V^*}^\perp =
\left(\begin{array}{cc}
0.125 & 0 \\
-0.114 & -0.172 \\
0.223 & -0.485 \\
0.833 & -0.102 \\
-0.131 & -0.253 \\
-0.339 & -0.462 \\
0.136 & 0.103 \\
0.218 & -0.0319 \\
-0.103 & 0.41 \\
-0.0721 & 0.0878 \\
0.102 & 0.406 \\
-0.0539 & 0.307 \\
\end{array}\right)
\end{align}
\begin{align} {\LARGE(12) \quad}
\mathcal L =
\left(\begin{array}{cccccccccc|}
-0.707 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0.406 & 0 & -0.414 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
-0.0811 & 0 & -0.159 & 0.981 & 0 & 0 & 0 & 0 & 0 & 0 \\
-0.406 & 0 & 0.414 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
-0.406 & 0 & -0.795 & -0.196 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & -1 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & -1 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & -1 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & -1 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -1 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \\
\end{array}\right)
\end{align}
\begin{align} {\LARGE(13) \quad}
T^T = \left(\begin{array}{cc|cccccccccc}
0.125 & 0 & -0.707 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
-0.114 & -0.172 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0.223 & -0.485 & 0.406 & 0 & -0.414 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0.833 & -0.102 & -0.0811 & 0 & -0.159 & 0.981 & 0 & 0 & 0 & 0 & 0 & 0 \\
-0.131 & -0.253 & -0.406 & 0 & 0.414 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
-0.339 & -0.462 & -0.406 & 0 & -0.795 & -0.196 & 0 & 0 & 0 & 0 & 0 & 0 \\
0.136 & 0.103 & 0 & 0 & 0 & 0 & -1 & 0 & 0 & 0 & 0 & 0 \\
0.218 & -0.0319 & 0 & 0 & 0 & 0 & 0 & -1 & 0 & 0 & 0 & 0 \\
-0.103 & 0.41 & 0 & 0 & 0 & 0 & 0 & 0 & -1 & 0 & 0 & 0 \\
-0.0721 & 0.0878 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -1 & 0 & 0 \\
0.102 & 0.406 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -1 & 0 \\
-0.0539 & 0.307 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \\
\end{array}\right)
\end{align}
tol = 1e-10;
prec = -log10(tol);
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} = round(T*A{i}/T, prec);
tB{i} = round(T*B{i}, prec);
tC{i} = round(C{i}/T, prec);
tD{i} = round(D{i}, prec);
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(14) \quad}
\bar A_{22}(\rho) =
\left(\begin{array}{cc|ccccccccc}
0.44\rho -2.2 & 0.074\rho +0.65 & 0.11\rho -0.45 & 0.099-0.21\rho & -0.21\rho -0.92 & -1.8\rho -10.0 & 1.2\rho -3.2 & 0.17\rho -3.5 & -3.2\rho -30.0 & 0.35-1.1\rho & -3.0\rho -16.0 \\
0.29\rho +1.5 & -0.25\rho -2.8 & 0.096\rho +1.3 & -0.24\rho -0.59 & -0.37\rho -10.0 & 0.085\rho +1.7 & 0.079\rho +2.5 & 0.077\rho -0.61 & -3.610^{-3}\rho -0.72 & 0.16\rho +1.8 & 0.18-0.067\rho \\ \hline
0.031\rho +0.43 & -0.11\rho -2.0 & 0.034\rho +0.15 & -0.013\rho -0.17 & -0.22\rho -3.1 & 0.037\rho +0.52 & 0.054\rho +0.74 & -0.013\rho -0.18 & -0.015\rho -0.21 & 0.038\rho +0.53 & 3.910^{-3}\rho +0.055 \\
0.38\rho +3.3 & 0.082\rho -0.2 & 0.12\rho +0.44 & -0.24\rho -2.5 & 1.0-0.39\rho & 0.1\rho +0.19 & 0.09\rho -0.35 & 0.11\rho +1.5 & 0.37-6.210^{-4}\rho & 0.2\rho +1.2 & -0.09\rho -1.1 \\
0.098\rho -0.13 & 0.01\rho +7.610^{-3} & 0.023\rho -0.017 & -0.074\rho -0.3 & -0.05\rho -0.039 & 0.018\rho -7.310^{-3} & 9.410^{-3}\rho +0.013 & 0.034\rho -0.057 & 4.510^{-3}\rho -0.014 & 0.046\rho -0.047 & 0.041-0.027\rho \\
0 & 0 & 0 & 0 & 0 & 0.61\rho +2.2 & -0.89\rho -0.19 & 1.9-0.21\rho & 1.1\rho +19.0 & 0.74\rho +1.2 & 1.4\rho +11.0 \\
0 & 0 & 0 & 0 & 0 & 2.1\rho +9.9 & 6.7-1.2\rho & 5.4-0.11\rho & 3.2\rho +28.0 & 1.1\rho -3.2 & 3.0\rho +12.0 \\
0 & 0 & 0 & 0 & 0 & -5.1\rho -21.0 & 3.9\rho -20.0 & 0.64\rho -16.0 & -7.2\rho -73.0 & 7.9-3.2\rho & -7.9\rho -29.0 \\
0 & 0 & 0 & 0 & 0 & 2.1\rho +8.9 & 6.6-1.6\rho & 5.8-0.26\rho & 2.9\rho +30.0 & 1.4\rho -2.7 & 3.1\rho +13.0 \\
0 & 0 & 0 & 0 & 0 & 5.6\rho +25.0 & 16.0-3.3\rho & 11.0-0.089\rho & 9.2\rho +81.0 & 3.2\rho -8.2 & 8.4\rho +38.0 \\
0 & 0 & 0 & 0 & 0 & -5.1\rho -24.0 & 3.9\rho -15.0 & 0.66\rho -13.0 & -7.1\rho -75.0 & 3.6-3.8\rho & -7.6\rho -36.0 \\
\end{array}\right)
\end{align}
\begin{align} {\LARGE(15) \quad}
\bar B(\rho) =
\left(\begin{array}{cc}
0.85 & 0.18 \\
0.29 & -1.5 \\ \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(16) \quad}
\bar C(\rho) = C = 10^{-3}
\left(\begin{array}{cc|cccccccccc}
1.5 & -0.05 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0.73 & -0.93 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
\end{array}\right)
\end{align}
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 = 1.0e+09 * -3.1174 -1.1826 -0.1368 -0.9343 -0.7637 -0.8061 -0.0000 -0.0000 -0.0000 -0.0000 ans = 1.0e+09 * -2.7092 -0.1686 -1.2891 -1.0916 -0.7694 -0.8056 -0.0000 -0.0000 -0.0000 -0.0000