Tartalomjegyzék

LPV output passivization

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.

1. lépés: Legyen egy LTI rendszermodell

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)
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
\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}

2. lépés: LPV modell

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 = B0*0;
Indices = [3];
B1(Indices) = rand(size(Indices));

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.385 & 0 & 0 & 0 & 0 & 0 \\ -0.084 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & -0.101 & 0 & 0.0327 & -0.18 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0.181 & 0 & 0.0121 & 0 \\ 0 & 0 & 0 & 0 & -0.156 & 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.724 & 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}

Observer design

Declare matrices

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

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*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*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
    ]

Solve the optimization problem

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

Validate solution

Zero dynamics of the obtained LPV

Construct LPV model matrices

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{:}];
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);

State transformation invariant feedback design

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 )

Transformed LPV system

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)

Check stability of the transformed system

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