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
];
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 ];
A_fh = @(rho) A0 + rho*A1;
B_fh = @(rho) B0 + rho*B1;
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.00
1 5.1e-01 2.2e-01 2.3e+00 1.00e+00 8.606546853e-01 0.000000000e+00 2.2e-01 0.01
2 1.7e-01 7.4e-02 1.4e+00 9.99e-01 2.946543699e-01 0.000000000e+00 7.4e-02 0.01
3 6.5e-02 2.8e-02 8.3e-01 9.96e-01 1.103446393e-01 0.000000000e+00 2.8e-02 0.01
4 2.2e-02 9.1e-03 4.7e-01 9.96e-01 3.656101357e-02 0.000000000e+00 9.1e-03 0.01
5 3.1e-03 1.3e-03 1.8e-01 9.91e-01 5.212144406e-03 0.000000000e+00 1.3e-03 0.01
6 1.3e-03 5.4e-04 1.1e-01 9.66e-01 2.208557457e-03 0.000000000e+00 5.4e-04 0.01
7 4.8e-04 2.0e-04 6.6e-02 9.39e-01 8.469565183e-04 0.000000000e+00 2.0e-04 0.02
8 1.2e-04 5.2e-05 3.0e-02 9.01e-01 2.254355678e-04 0.000000000e+00 5.1e-05 0.02
9 3.7e-05 1.6e-05 1.3e-02 7.30e-01 7.582956194e-05 0.000000000e+00 1.6e-05 0.02
10 1.0e-05 4.2e-06 3.0e-03 1.89e-01 1.920705892e-05 0.000000000e+00 4.2e-06 0.02
11 2.4e-06 1.0e-06 3.8e-04 -5.94e-01 -1.077878747e-04 0.000000000e+00 1.0e-06 0.02
12 6.0e-07 2.5e-07 4.7e-05 -9.13e-01 -6.511558896e-04 0.000000000e+00 2.5e-07 0.02
13 1.9e-07 8.1e-08 9.2e-06 -9.16e-01 -1.857025229e-03 0.000000000e+00 8.1e-08 0.02
14 4.8e-08 2.0e-08 1.1e-06 -1.05e+00 -8.810702944e-03 0.000000000e+00 2.0e-08 0.02
15 1.2e-08 5.2e-09 1.5e-07 -9.39e-01 -2.916081443e-02 0.000000000e+00 5.2e-09 0.02
16 3.2e-09 1.4e-09 2.0e-08 -9.59e-01 -1.122226036e-01 0.000000000e+00 1.4e-09 0.03
17 9.0e-10 3.8e-10 3.1e-09 -9.48e-01 -3.752172278e-01 0.000000000e+00 3.8e-10 0.03
18 2.4e-10 1.0e-10 4.4e-10 -9.59e-01 -1.367181983e+00 0.000000000e+00 1.0e-10 0.03
19 5.5e-11 2.3e-11 5.1e-11 -9.16e-01 -5.206768002e+00 0.000000000e+00 2.3e-11 0.03
20 1.2e-11 5.0e-12 5.2e-12 -9.16e-01 -2.238294604e+01 0.000000000e+00 5.0e-12 0.03
21 2.6e-12 1.1e-12 5.4e-13 -9.58e-01 -1.018078662e+02 0.000000000e+00 1.1e-12 0.03
22 6.1e-13 2.8e-13 6.6e-14 -9.33e-01 -3.848821263e+02 0.000000000e+00 2.6e-13 0.03
23 1.3e-13 8.9e-14 6.9e-15 -9.10e-01 -1.585165169e+03 0.000000000e+00 5.5e-14 0.04
24 2.2e-14 6.7e-13 4.9e-16 -9.07e-01 -8.706179538e+03 0.000000000e+00 9.1e-15 0.04
25 2.1e-15 4.9e-13 9.7e-18 -9.56e-01 -1.131541259e+05 0.000000000e+00 6.5e-16 0.04
26 1.4e-15 3.8e-13 6.4e-18 -9.89e-01 -1.496285797e+05 0.000000000e+00 4.9e-16 0.04
27 1.0e-15 3.1e-13 2.3e-18 -9.92e-01 -2.919815748e+05 0.000000000e+00 2.5e-16 0.05
28 1.0e-15 3.1e-13 2.3e-18 -9.96e-01 -2.919815748e+05 0.000000000e+00 2.5e-16 0.05
Optimizer terminated. Time: 0.06
MOSEK DUAL INFEASIBILITY REPORT.
Problem status: The problem is dual infeasible
Interior-point solution summary
Problem status : DUAL_INFEASIBLE
Solution status : DUAL_INFEASIBLE_CER
Primal. obj: -8.6357352583e-06 nrm: 1e+01 Viol. con: 2e-13 barvar: 2e-16
Optimizer summary
Optimizer - time: 0.06
Interior-point - iterations : 29 time: 0.05
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.2549
solvertime: 0.0641
info: 'Infeasible problem (MOSEK)'
problem: 1
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
| ID| Constraint| Primal residual| Dual residual|
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
| #1| Matrix inequality| -1.5383e-16| 1.542e-16|
| #2| Matrix inequality| -4.9759e-16| 1.8597e-16|
| #3| Matrix inequality| -0.0001| 8.9219e-17|
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
[ OK ] M2(0) <= 0 [ OK ] M2(0.1) <= 0 [ OK ] M2(-0.1) <= 0 [ OK ] P > 0 [ OK ] A - LC < 0 [FAILED] norm(D1) > 1e-3 2.8713e-11
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{:}];
syms rho real
if ~exist('pcz_latex_output','file')
pcz_latex_output = @(varargin) [];
end
C_mtp = round(log10(norm(C{1})));
pcz_latex_output(vpa(A_fh(rho),2), 'disp_mode', 2, 'label', '\nA(\\rho) = ');
pcz_latex_output(vpa(B_fh(rho),2), 'disp_mode', 2, 'label', '\nB(\\rho) = ');
pcz_latex_output(vpa(C_fh(rho) / 10^C_mtp,2), 'disp_mode', 2, 'label', ['\nC(\\rho) = C = 10^{' num2str(C_mtp) '}' ]);
Nálam az $E_c$ most az $Im(B)$.
Im_B = IMA([B{1} B{2}]);
pcz_num2str_latex_output(Im_B, 'label', '\mathrm{Im}\big(B(\rho)\big) = ', 'esc', 0);
Ker_C = INTS(KER(C{1}), KER(C{2}));
pcz_num2str_latex_output(Ker_C, 'label', '\\mathrm{Ker}(C) = ');
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\}$.
pcz_num2str_latex_output(V, 'label', '\mathcal V^* = ', 'esc', 0);
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) ]';
pcz_num2str_latex_output(ORTCO(V), 'label', [ '{\mathcal V^*}^\perp = ' newline ' ' ], 'esc', 0)
pcz_num2str_latex_output(L, 'label', [ '\mathcal L = ' newline ' ' ], 'vline', dim_V, 'esc', 0)
pcz_num2str_latex_output(T', 'label', 'T^T = ', 'vline', dim_Ker_V )
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;
syms rho real
A_rho = tA_fh(rho);
A22 = A_rho(dim_Ker_V:end,dim_Ker_V:end);
pcz_latex_output(vpa(A22,2), 'disp_mode', 2, 'label', '\n\\bar A_{22}(\\rho) = ', 'vline', dim_Ker_V, 'hline', dim_Ker_V)
pcz_latex_output(vpa(tB_fh(rho),2), 'disp_mode', 2, 'label', '\n\\bar B(\\rho) = ', 'hline', dim_Ker_V)
pcz_latex_output(vpa(tC_fh(rho) / 10^C_mtp,2), 'disp_mode', 2, 'label', ['\n\\bar C(\\rho) = C = 10^{' num2str(C_mtp) '}'], 'vline', dim_Ker_V)
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