Matlab Exercise

Chapter 2: システムの状態空間表現

---

2.2 線形システムの状態空間表現

2.2.2 同値変換(equivalence transformation)

$\left\{\begin{array}{l}\dot{\overline{x}}\left(t\right)=\overline{A}\overline{x}\left(t\right)+\overline{B}u\left(t\right)\\ y\left(t\right)=\overline{C}\overline{x}\left(t\right)+Du\left(t\right)\end{array}\right., \overline{x}\left(t\right)=Tx\left(t\right)\in {\mathrm{ℝ}}^n, \overline{A}=TA{T}^{-1}\in {\mathrm{ℝ}}^{n\times n},$

$\overline{B}=TB\in {\mathrm{ℝ}}^{n\times p}, \overline{C}=C{T}^{-1}\in {\mathrm{ℝ}}^{p\times n} (2.1)$

$T$ is regular matrix (invertible matrix, 可逆矩阵), any regular matrix can be applied in Equation (2.1) , thus, nember of state-spaces is countless as transform matrix $T$ changes.

例2.5 state-space of RCL circuit and its equivalence transformation

example2_5.m

%% vc as states
syms vc(t) i(t) C x vl(t) R vr(t) L vin u
eq1 = vin == vl(t)+vr(t)+vc(t);
eq2 = vl(t) == L*diff(i(t),t);
eq3 = vr(t) == R*i(t);
eq4 = diff(vc(t),t) == diff(1/C*int(i(t),t), t);
eq4x = subs(eq4,i,x);
cancelI = solve(eq4x,x);
subi2 = subs(eq2,i,cancelI);
subi3 = subs(eq3,i,cancelI);
% rhs is available after version 2017a
% subi2 = rhs(subs(eq2,i,cancelI));
% subi3 = rhs(subs(eq3,i,cancelI));
eq = subs(subs(eq1,vl,C*L*diff(vc(t), t, t)),vr, C*R*diff(vc(t), t));
[vf,Yvar] = odeToVectorField(eq);
F_vf = matlabFunction(vf);
x_vf = F_vf(C,L,R,Yvar,u);
mat_A = equationsToMatrix(x_vf ==0, Yvar);
mat_B = equationsToMatrix(x_vf ==0, u);

%% charge as states
syms q(t)
eq4q = vc(t) == 1/C*int(i(t),t);
eq5 = diff(q(t),t) == diff(int(i(t),t),t);
eq5x = subs(eq5,i,x);
canceliq = solve(eq5x,x);
subi2q = subs(eq2,i,canceliq);
subi3q = subs(eq3,i,canceliq);
subi4q = subs(eq4q,i,canceliq);
eqq = subs(subs(subs(eq1,vl,L*diff(q(t), t, t)),vr, R*diff(q(t), t)),vc,q(t)/C);
[vfq,Yvarq] = odeToVectorField(eqq);
F_vfq = matlabFunction(vfq);
x_vfq = F_vfq(C,L,R,Yvarq,u);
mat_Aq = equationsToMatrix(x_vfq ==0, Yvarq);
mat_Bq = equationsToMatrix(x_vfq ==0, u);
fout = subs(sym('y=1/C*int(i(t),t)'),i,canceliq);

問題2.2 equivalence transformation

problem2_2.m
ans: (1)

syms x1(t) x2(t) x3(t) x4(t) u;
syms z1(t) z2(t);
eq1 = 0.5*diff(z1(t),t,t) == u-2*(z1(t)-z2(t))-(diff(z1(t),t)-diff(z2(t),t));
eq2 = diff(z2(t),t,t) == 2*(z1(t)-z2(t))+(diff(z1(t),t)-diff(z2(t),t));

[vf,Yvar] = odeToVectorField([eq1,eq2]);

rearr = matlabFunction(Yvar);
z_state_rearr = rearr(diff(z1,t),diff(z2,t),z1,z2);%rearr=@(Dz1,Dz2,z1,z2)...
z_state = [z1;diff(z1,t);z2;diff(z2,t)];
[i,j] = find(z_state*transpose(1./z_state_rearr)==1);
TrM = zeros(4);
TrM(sub2ind([4,4],i,j)) = 1;
sort_vf = TrM*vf;
F_vf = matlabFunction(sort_vf);
z_vf = F_vf(Yvar,u);
mat_A = double(equationsToMatrix(z_vf ==0, TrM*Yvar));
mat_B = double(equationsToMatrix(z_vf ==0, u));

syms x11 x22 x33 x44 z11 z22 dz11 dz22
dz_state = [z11; dz11; z22; dz22];
fun = formula(mat_A*dz_state+mat_B*u);

dz1 = fun(1);
ddz1 = fun(2);
dz2 = fun(3);
ddz2 = fun(4);
dx1 = dz1;
dx2 = ddz1;
dx3 = dz1-dz2;
dx4 = ddz1-ddz2;

funxz = [x11 == z11, x22 == dz11, x33 == z11 - z22, x44 == dz11 - dz22];
[cancelz11,canceldz11,cancelz22,canceldz22]=solve(funxz,[z11,dz11,z22,dz22]);
rhsfunx = subs([dx1,dx2,dx3,dx4],[z11,dz11,z22,dz22],[cancelz11,canceldz11,cancelz22,canceldz22]);
rhsfunxt = subs(rhsfunx,[x11,x22,x33,x44],[x1,x2,x3,x4]);
odext = formula([diff(x1,t),diff(x2,t),diff(x3,t),diff(x4,t)])==rhsfunxt;
[vfx,Yvarx] = odeToVectorField(odext);

rearrx = matlabFunction(Yvarx);
x_state_rearr = rearrx(x1,x2,x3,x4);
x_state = [x1;x2;x3;x4];
[ix,jx] = find(x_state*transpose(1./x_state_rearr)==1);
TrMx = zeros(4);
TrMx(sub2ind([4,4],ix,jx)) = 1;
sort_vfx = TrMx*vfx;
F_vfx = matlabFunction(sort_vfx);
x_vf = F_vfx(Yvarx,u);
mat_Ax = double(equationsToMatrix(x_vf ==0, TrMx*Yvarx));
mat_Bx = double(equationsToMatrix(x_vf ==0, u));

Appendix A. 矩阵的等价(equivalence)，相似(similarity)，合同(conguence)

Definition:
等价矩阵： 有相同的秩。
相似矩阵： 有相同的特征值和秩， $P^{-1}AP=B$
合同矩阵： 有相同的正负惯性指数，即正负特征值个数相同， $P^{T}AP=B$

与相似矩阵相关的概念：特征值，矩阵对角化，约旦标准型，齐次线性方程。
与合同矩阵相关的概念：二次型，标准型。