Teljes Matlab script kiegészítő függvényekkel.
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.
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,B,C,D] = deal(sys.a, sys.b, sys.c, sys.d);
tol = 1e-10;
prec = -log10(tol);
A0 = round(A0,prec);
B = round(B,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} $$
A1 = A0;
A1(abs(A0) < 1) = 0;
A1 = A1/100;
B0 = B;
B1 = B*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.0537 & 0 & 0 & 0 & 0 & 0 \\
0.0187 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & -0.0135 & 0 & -0.0155 & -0.0155 \\
0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & -0.0319 & 0 & -0.0181 & 0 \\
0 & 0 & 0 & 0 & 0.019 & 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,B,linspace(-5,-6,n_x));
C1 = sdpvar(n_yp,n_x,'full');
D1 = sdpvar(n_yp,n_y,'full');
$$ \begin{align} &\hat A(\rho) = \pmqty{ A(\rho) & 0 \\ 0 & A(\rho) - L C },~ \hat B = \pmqty{B \\ 0} \\ &\hat 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*C
];
Bo = @(rho) [B_fh(rho) ; zeros(n_x,n_u)];
Co = [D1*C+C1 -C1];
Do = zeros(n_yp,n_r);
$$ \begin{aligned} \hat A_c(\rho) = \begin{pmatrix} A(\rho)-B(\rho) K & B(\rho) K \\ 0 & A(\rho) - L C \end{pmatrix},~ \hat 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*C
];
Bc = @(rho) Bo(rho)*G;
W = eye(n_yp);
Q = sdpvar(n_x,n_x,'symmetric');
S = sdpvar(n_x,n_x,'symmetric');
N = sdpvar(n_x,n_y,'full');
P = blkdiag(Q,S);
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*C-C'*N'
];
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(-1) <= 0
M2(1) <= 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)
Problem Name : Objective sense : min Type : CONIC (conic optimization problem) Constraints : 70 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 : 70 Cones : 0 Scalar variables : 0 Matrix variables : 3 Integer variables : 0 Optimizer - threads : 4 Optimizer - solved problem : the primal Optimizer - Constraints : 70 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 : 2485 after factor : 2485 Factor - dense dim. : 0 flops : 4.63e+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.1e-01 2.2e-01 2.3e+00 1.00e+00 8.673126699e-01 0.000000000e+00 2.2e-01 0.02 2 1.7e-01 7.3e-02 1.3e+00 9.99e-01 2.921680370e-01 0.000000000e+00 7.3e-02 0.02 3 6.5e-02 2.8e-02 8.3e-01 9.96e-01 1.107970153e-01 0.000000000e+00 2.8e-02 0.02 4 2.4e-02 1.0e-02 5.0e-01 9.96e-01 4.030266395e-02 0.000000000e+00 1.0e-02 0.02 5 4.9e-03 2.1e-03 2.3e-01 9.93e-01 8.390005047e-03 0.000000000e+00 2.1e-03 0.02 6 9.6e-04 4.0e-04 9.7e-02 9.80e-01 1.639123284e-03 0.000000000e+00 4.0e-04 0.03 7 3.8e-04 1.6e-04 6.0e-02 9.66e-01 6.547062625e-04 0.000000000e+00 1.6e-04 0.03 8 1.0e-04 4.3e-05 3.1e-02 9.73e-01 1.788249199e-04 0.000000000e+00 4.3e-05 0.03 9 3.4e-05 1.4e-05 1.7e-02 9.63e-01 6.007841566e-05 0.000000000e+00 1.4e-05 0.04 10 1.1e-05 4.5e-06 9.3e-03 9.36e-01 1.968189299e-05 0.000000000e+00 4.5e-06 0.04 11 3.6e-06 1.5e-06 5.0e-03 8.79e-01 6.959485711e-06 0.000000000e+00 1.5e-06 0.04 12 1.1e-06 4.6e-07 2.3e-03 7.87e-01 2.279750611e-06 0.000000000e+00 4.6e-07 0.04 13 4.0e-07 1.7e-07 1.1e-03 6.12e-01 9.464609679e-07 0.000000000e+00 1.7e-07 0.04 14 1.1e-07 4.8e-08 4.0e-04 4.51e-01 2.573487329e-07 0.000000000e+00 4.8e-08 0.05 15 4.2e-08 1.8e-08 1.6e-04 2.56e-01 5.300163191e-08 0.000000000e+00 1.8e-08 0.05 16 1.2e-08 4.9e-09 4.9e-05 1.78e-01 -8.163789148e-08 0.000000000e+00 4.9e-09 0.05 17 4.2e-09 1.8e-09 1.8e-05 8.39e-02 -1.349185641e-07 0.000000000e+00 1.8e-09 0.05 18 1.1e-09 4.8e-10 5.2e-06 7.17e-02 -1.643388631e-07 0.000000000e+00 4.8e-10 0.05 19 4.2e-10 1.8e-10 1.9e-06 1.58e-02 -1.876681766e-07 0.000000000e+00 1.8e-10 0.06 20 1.1e-10 4.7e-11 5.3e-07 3.98e-02 -1.874660894e-07 0.000000000e+00 4.7e-11 0.06 21 4.3e-11 1.8e-11 2.0e-07 -6.86e-03 -2.041125985e-07 0.000000000e+00 1.8e-11 0.06 22 1.1e-11 4.8e-12 5.4e-08 2.71e-02 -1.960984933e-07 0.000000000e+00 4.8e-12 0.07 23 4.4e-12 1.9e-12 2.0e-08 -9.23e-04 -2.051474785e-07 0.000000000e+00 1.9e-12 0.07 24 1.2e-12 5.1e-13 5.6e-09 1.93e-02 -1.992451236e-07 0.000000000e+00 5.0e-13 0.07 25 4.2e-13 1.8e-13 2.0e-09 7.88e-04 -2.045000820e-07 0.000000000e+00 1.8e-13 0.08 26 1.2e-13 5.5e-14 5.9e-10 2.44e-02 -1.961453373e-07 0.000000000e+00 5.2e-14 0.08 27 4.0e-14 2.2e-14 1.8e-10 -1.63e-02 -2.094143590e-07 0.000000000e+00 1.7e-14 0.08 28 1.2e-14 6.2e-15 5.8e-11 5.45e-02 -1.846388577e-07 0.000000000e+00 5.0e-15 0.08 29 6.7e-15 1.1e-14 2.0e-11 -4.15e-02 -2.150437599e-07 0.000000000e+00 1.8e-15 0.08 30 1.2e-14 4.2e-15 4.6e-12 1.55e-02 -1.967879638e-07 0.000000000e+00 4.2e-16 0.08 31 1.3e-14 8.9e-15 2.0e-12 1.30e-02 -2.058178671e-07 0.000000000e+00 1.8e-16 0.08 32 5.5e-15 7.7e-15 5.3e-13 -2.34e-02 -2.056777324e-07 0.000000000e+00 4.9e-17 0.08 33 7.5e-15 4.3e-15 1.9e-13 1.39e-02 -2.005825246e-07 0.000000000e+00 1.7e-17 0.08 34 1.1e-14 1.1e-14 5.9e-14 1.31e-02 -2.016908340e-07 0.000000000e+00 5.3e-18 0.09 35 5.8e-15 7.3e-15 1.7e-14 -2.07e-02 -2.005799148e-07 0.000000000e+00 1.5e-18 0.09 36 1.5e-14 5.1e-15 7.2e-15 3.45e-02 -1.982411837e-07 0.000000000e+00 6.5e-19 0.09 37 1.8e-14 2.9e-15 1.6e-15 -1.38e-02 -2.052754923e-07 0.000000000e+00 1.5e-19 0.09 38 1.2e-14 3.9e-15 8.1e-16 4.01e-02 -1.987062507e-07 0.000000000e+00 7.3e-20 0.09 39 3.5e-15 5.6e-15 2.4e-16 6.90e-03 -2.059511331e-07 0.000000000e+00 2.2e-20 0.09 40 1.0e-15 5.6e-15 6.8e-17 -1.93e-02 -2.091418935e-07 0.000000000e+00 6.4e-21 0.10 41 3.6e-16 1.3e-14 2.5e-17 2.12e-02 -1.999897939e-07 0.000000000e+00 2.3e-21 0.10 42 1.2e-16 4.0e-15 8.0e-18 2.57e-02 -2.038122505e-07 0.000000000e+00 7.4e-22 0.10 43 4.1e-17 4.9e-15 2.9e-18 2.08e-02 -1.964082699e-07 0.000000000e+00 2.6e-22 0.10 44 1.6e-17 6.0e-15 1.1e-18 4.14e-02 -1.997433769e-07 0.000000000e+00 9.8e-23 0.10 45 4.2e-18 8.0e-15 2.9e-19 -5.26e-03 -2.054041850e-07 0.000000000e+00 2.7e-23 0.11 46 1.3e-18 3.1e-15 9.2e-20 6.94e-03 -2.031176491e-07 0.000000000e+00 8.4e-24 0.11 47 4.5e-19 7.1e-15 3.1e-20 -1.47e-02 -2.081756359e-07 0.000000000e+00 2.8e-24 0.11 48 6.4e-18 5.8e-15 9.9e-21 9.84e-03 -2.004502391e-07 0.000000000e+00 9.0e-25 0.11 49 3.1e-17 8.0e-15 8.3e-21 -2.73e-03 -2.008871074e-07 0.000000000e+00 7.5e-25 0.11 50 3.7e-17 8.1e-15 8.1e-21 -6.46e-04 -2.009007380e-07 0.000000000e+00 7.3e-25 0.12 51 3.8e-17 7.5e-15 8.1e-21 4.31e-04 -2.009016196e-07 0.000000000e+00 7.3e-25 0.12 52 3.8e-17 7.5e-15 8.1e-21 1.34e-03 -2.009016196e-07 0.000000000e+00 7.3e-25 0.12 Optimizer terminated. Time: 0.13 Interior-point solution summary Problem status : PRIMAL_AND_DUAL_FEASIBLE Solution status : OPTIMAL Primal. obj: -2.0090161957e-07 nrm: 2e+09 Viol. con: 3e-07 barvar: 2e-07 Dual. obj: 0.0000000000e+00 nrm: 3e+09 Viol. con: 0e+00 barvar: 4e-06 Optimizer summary Optimizer - time: 0.13 Interior-point - iterations : 53 time: 0.13 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.2536 solvertime: 0.1449 info: 'Successfully solved (MOSEK)' problem: 0 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ | ID| Constraint| Primal residual| Dual residual| ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ | #1| Matrix inequality| -1.4044e-07| 5.5714e-17| | #2| Matrix inequality| -5.225e-08| 8.5553e-17| | #3| Matrix inequality| 2.4179e-05| 5.3487e-17| ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Q = value(Q);
S = value(S);
N = value(N);
L = S\N;
P = value(P);
Co = value(Co);
C1 = value(C1);
D1 = value(D1);
% value(AcP_PAc2(0.1) - AcP_PAc(0.1))
Bc0 = Bc(0);
Bo0 = Bc(0);
Ac0 = Ac(L,0);
Ao0 = Ao(L,0);
M2_val = value(M2(0));
if pcz_info(all(real(eig(M2_val)) <= 1e-5), 'M2(0) <= 0')
eig(M2_val)'
fprintf('larges eigen value: %f\n\n', max(eig(M2_val)))
end
M2_val = value(M2(0.1));
if pcz_info(all(real(eig(M2_val)) <= 1e-5), 'M2(0.1) <= 0')
eig(M2_val)'
fprintf('larges eigen value: %f\n\n', max(eig(M2_val)))
end
M2_val = value(M2(-0.1));
if pcz_info(all(real(eig(M2_val)) <= 1e-5), 'M2(-0.1) <= 0')
eig(M2_val)'
fprintf('larges eigen value: %f\n\n', max(eig(M2_val)))
end
if pcz_info(all(real(eig(P)) > 0), 'P > 0')
eig(P)'
end
if pcz_info(all(real(eig(A0 - L*C)) < 0), 'A - LC < 0')
eig(A0 - L*C)
end
if pcz_info(norm(D1) > 1e-3, 'norm(D1) > 1e-3')
disp(norm(D1))
end
% The resulting closed-loop system
sys_CLS = minreal(ss(Ac0,Bc0,Co,Do),[],0);
sys_OLS = minreal(ss(Ao0,Bo0,Co,Do),[],0);
sys_CLS = minreal(tf(sys_CLS),[],0);
sys_OLS = minreal(tf(sys_OLS),[],0);
if pcz_info(all(real(tzero(sys_OLS)) < 0), 'Zeros of the open loop system are stable')
tzero(sys_OLS)
end
if pcz_info(all(real(tzero(sys_CLS)) < 0), 'Zeros of the closed loop system are stable')
tzero(sys_CLS)
end
[POLES,ZEROS] = pzmap(sys_OLS);
[POLES,ZEROS] = pzmap(sys_CLS);
[ 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,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.054\rho -5.4 & -0.93 & -0.46 & 0.093 & 0.46 & 0.46 & 0 & 0 & 0 & 0 & 0 & 0 \\
0.019\rho +1.9 & -0.32 & -0.19 & 0.038 & 0.19 & 0.19 & 0 & 0 & 0 & 0 & 0 & 0 \\
-0.29 & -0.73 & -0.014\rho -1.4 & 0.29 & -0.015\rho -1.5 & -0.015\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.032\rho -3.2 & 0.64 & -0.018\rho -1.8 & 0.19 & 0 & 0 & 0 & 0 & 0 & 0 \\
0.29 & 0.73 & 0.1 & -0.02 & 0.019\rho +1.9 & -0.1 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 1.8-0.054\rho & -0.61 & 1.6 & 18.0 & 1.4 & -10.0 \\
0 & 0 & 0 & 0 & 0 & 0 & 0.019\rho +10.0 & -1.3 & 1.5 & 19.0 & 2.1 & -12.0 \\
0 & 0 & 0 & 0 & 0 & 0 & -21.0 & 0.72 & -0.014\rho -5.9 & -47.0 & -0.015\rho -5.7 & 29.0-0.015\rho \\
0 & 0 & 0 & 0 & 0 & 0 & 8.1 & -0.29 & 2.1 & 19.0 & 1.8 & -12.0 \\
0 & 0 & 0 & 0 & 0 & 0 & 26.0 & -1.4 & 2.3-0.032\rho & 59.0 & 3.6-0.018\rho & -38.0 \\
0 & 0 & 0 & 0 & 0 & 0 & 23.0 & -1.6 & 4.8 & 52.0 & 0.019\rho +7.0 & -34.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.23 & -0.2 & 0.34 & 1.4 & -0.13 & -0.6 & 0.23 & 0.34 & -0.15 & -0.15 & 0.11 & -0.088 \\
0.095 & 0.065 & 0.68 & 0.82 & 0.095 & 0.22 & 6.710^{-4} & 0.26 & -0.49 & -0.17 & -0.3 & -0.37 \\
\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.991 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
-0.0193 & 0.979 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0.0173 & -0.0591 & -0.84 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0.113 & 0.0889 & 0.254 & 0.44 & 0 & 0 & 0 & 0 & 0 & 0 \\
-0.0137 & -0.0363 & 0.0783 & 0.0822 & 0.977 & 0 & 0 & 0 & 0 & 0 \\
-0.0585 & -0.131 & 0.232 & 0.482 & -0.178 & 0.515 & 0 & 0 & 0 & 0 \\
0.0206 & 0.0359 & -0.0374 & -0.228 & 0.0532 & 0.468 & -0.832 & 0 & 0 & 0 \\
0.0253 & 0.0099 & 0.1 & -0.468 & 0.0337 & 0.459 & 0.405 & 0.588 & 0 & 0 \\
-0.0041 & 0.0573 & -0.269 & 0.404 & 0.0489 & 0.124 & 0.0178 & 0.347 & -0.669 & 0 \\
-0.0103 & 0.004 & -0.0757 & 0.237 & -0.0063 & -0.161 & -0.154 & 0.477 & 0.455 & -0.661 \\
0.0159 & 0.0667 & -0.196 & 0.0418 & 0.077 & 0.492 & 0.341 & -0.528 & 0.145 & -0.405 \\
-0.0006 & 0.047 & -0.205 & 0.275 & 0.0428 & 0.148 & 0.0571 & 0.162 & 0.57 & 0.631 \\
\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.991 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ -0.0193 & 0.979 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0.0173 & -0.0591 & -0.84 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0.113 & 0.0889 & 0.254 & 0.44 & 0 & 0 & 0 & 0 & 0 & 0 \\ -0.0137 & -0.0363 & 0.0783 & 0.0822 & 0.977 & 0 & 0 & 0 & 0 & 0 \\ -0.0585 & -0.131 & 0.232 & 0.482 & -0.178 & 0.515 & 0 & 0 & 0 & 0 \\ 0.0206 & 0.0359 & -0.0374 & -0.228 & 0.0532 & 0.468 & -0.832 & 0 & 0 & 0 \\ 0.0253 & 0.0099 & 0.1 & -0.468 & 0.0337 & 0.459 & 0.405 & 0.588 & 0 & 0 \\ -0.0041 & 0.0573 & -0.269 & 0.404 & 0.0489 & 0.124 & 0.0178 & 0.347 & -0.669 & 0 \\ -0.0103 & 0.004 & -0.0757 & 0.237 & -0.0063 & -0.161 & -0.154 & 0.477 & 0.455 & -0.661 \\ 0.0159 & 0.0667 & -0.196 & 0.0418 & 0.077 & 0.492 & 0.341 & -0.528 & 0.145 & -0.405 \\ -0.0006 & 0.047 & -0.205 & 0.275 & 0.0428 & 0.148 & 0.0571 & 0.162 & 0.57 & 0.631 \\ \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.136 & 0 \\
-0.141 & -0.145 \\
0.126 & -0.524 \\
0.822 & -0.214 \\
-0.1 & -0.146 \\
-0.428 & -0.462 \\
0.151 & 0.0929 \\
0.185 & -0.117 \\
-0.0301 & 0.417 \\
-0.0756 & 0.102 \\
0.116 & 0.336 \\
-0.0046 & 0.322 \\
\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.136 & 0 & -0.707 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
-0.141 & -0.145 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0.126 & -0.524 & 0.406 & 0 & -0.414 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0.822 & -0.214 & -0.0811 & 0 & -0.159 & 0.981 & 0 & 0 & 0 & 0 & 0 & 0 \\
-0.1 & -0.146 & -0.406 & 0 & 0.414 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
-0.428 & -0.462 & -0.406 & 0 & -0.795 & -0.196 & 0 & 0 & 0 & 0 & 0 & 0 \\
0.151 & 0.0929 & 0 & 0 & 0 & 0 & -1 & 0 & 0 & 0 & 0 & 0 \\
0.185 & -0.117 & 0 & 0 & 0 & 0 & 0 & -1 & 0 & 0 & 0 & 0 \\
-0.0301 & 0.417 & 0 & 0 & 0 & 0 & 0 & 0 & -1 & 0 & 0 & 0 \\
-0.0756 & 0.102 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -1 & 0 & 0 \\
0.116 & 0.336 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -1 & 0 \\
-0.0046 & 0.322 & 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(\rho) =
\left(\begin{array}{cc|cccccccccc}
-9.410^{-3}\rho -1.8 & 9.710^{-3}\rho +0.39 & 0.012\rho +1.3 & 8.410^{-5}\rho -0.64 & -6.810^{-3}\rho -0.43 & 9.910^{-3}\rho +1.7 & 4.110^{-3}\rho -5.3 & 0.11-2.910^{-3}\rho & 7.610^{-3}\rho -0.56 & 1.710^{-3}\rho -13.0 & 3.910^{-3}\rho -1.1 & -2.710^{-3}\rho -8.2 \\
-7.410^{-4}\rho -0.016 & -0.021\rho -1.8 & 3.710^{-3}\rho +0.43 & -3.110^{-3}\rho -0.31 & 1.010^{-3}\rho +0.054 & -3.410^{-3}\rho -0.6 & 5.210^{-3}\rho -7.4 & 2.310^{-3}\rho +0.84 & 7.810^{-3}\rho -0.81 & -2.010^{-3}\rho -18.0 & -5.510^{-4}\rho -1.8 & 1.410^{-4}\rho -12.0 \\ \hline
0.12\rho +12.0 & 3.410^{-3}\rho +0.48 & -0.031\rho -2.5 & 0.017\rho +1.6 & -0.016\rho -1.7 & -0.1\rho -11.0 & 0.018\rho +1.9 & 0.022\rho +2.2 & -2.210^{-3}\rho -0.17 & -8.610^{-3}\rho -0.89 & 0.015\rho +1.6 & -5.310^{-4}\rho -0.097 \\
0.045\rho +3.7 & 1.710^{-3}\rho +0.14 & -0.018\rho -1.9 & 6.510^{-3}\rho +0.22 & -6.110^{-3}\rho -0.51 & -0.04\rho -3.3 & 6.910^{-3}\rho +0.57 & 8.010^{-3}\rho +0.66 & -6.210^{-4}\rho -0.052 & -3.210^{-3}\rho -0.26 & 5.810^{-3}\rho +0.48 & -3.510^{-4}\rho -0.029 \\
-3.210^{-3}\rho -0.65 & 0.035\rho +3.8 & -6.010^{-4}\rho -0.13 & 4.710^{-3}\rho +0.45 & -0.028\rho -3.0 & 5.810^{-3}\rho +1.0 & 2.810^{-3}\rho +0.25 & -4.710^{-3}\rho -0.56 & 0.015\rho +1.6 & 3.910^{-3}\rho +0.43 & 0.011\rho +1.2 & -0.011\rho -1.2 \\
4.110^{-3}\rho +0.025 & 5.710^{-3}\rho -0.15 & 7.910^{-4}\rho +4.810^{-3} & 1.410^{-3}\rho -0.017 & -5.210^{-3}\rho -0.29 & -3.010^{-3}\rho -0.04 & 1.210^{-3}\rho -9.710^{-3} & 9.810^{-5}\rho +0.022 & 2.310^{-3}\rho -0.061 & 2.710^{-4}\rho -0.017 & 2.410^{-3}\rho -0.046 & 0.047-1.810^{-3}\rho \\
0 & 0 & 0 & 0 & 0 & 0 & 1.8-0.054\rho & -0.61 & 1.6 & 18.0 & 1.4 & 10.0 \\
0 & 0 & 0 & 0 & 0 & 0 & 0.019\rho +10.0 & -1.3 & 1.5 & 19.0 & 2.1 & 12.0 \\
0 & 0 & 0 & 0 & 0 & 0 & -21.0 & 0.72 & -0.014\rho -5.9 & -47.0 & -0.015\rho -5.7 & 0.015\rho -29.0 \\
0 & 0 & 0 & 0 & 0 & 0 & 8.1 & -0.29 & 2.1 & 19.0 & 1.8 & 12.0 \\
0 & 0 & 0 & 0 & 0 & 0 & 26.0 & -1.4 & 2.3-0.032\rho & 59.0 & 3.6-0.018\rho & 38.0 \\
0 & 0 & 0 & 0 & 0 & 0 & -23.0 & 1.6 & -4.8 & -52.0 & -0.019\rho -7.0 & -34.0 \\
\end{array}\right)
\end{align}
\begin{align} {\LARGE(15) \quad}
\bar B(\rho) =
\left(\begin{array}{cc}
0.83 & 0.052 \\
0.14 & -1.3 \\ \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.7 & -0.25 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0.7 & -1.1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
\end{array}\right)
\end{align}