Matlab Exercise

Chapter 1: 古典制御理論から現代制御理論へ

---

→ １．１高次システムに対する古典制御理論の限界

例１．１ transfer function representation

example1_1.m

syms s y u z
eq1 = 0.5*s^2*z==u-2*(z-y)-(s*z-s*y);
eq2 = s^2*y==2*(z-y)+(s*z-s*y);         %input simultaneous equations

cancelZ = solve(eq2, z);                %cancel Z from eq2
uy = subs(eq1, z, cancelZ);             %substitute Z into eq1
Ps = subs(solve(uy, y),u,1);            %symbolic transfer function
[symnum,symden] = numden(Ps);           %extract num den
num = sym2poly(symnum);                 
den = sym2poly(symden);                 %matrix format

tf = tf(num,den);                       %transfer function
%num = cell2mat(tf.numerator(1,1));     
%den = cell2mat(tf.denominator(1,1));   %cell to matrix

example1_1.m shows how to generate transfer function from frequency domain equations (1.1)

$\left\{\begin{array}{l}0.5s^{2}y\left(s\right)=2\left(z_1\right(s)-y(s\left)\right)+\left(s{z}_1\right(s)-sy(s\left)\right)\\ s^2{z}_1\left(s\right)=u\left(s\right)-2\left(z_1\right(s)-y(s\left)\right)-\left(s{z}_1\right(s)-sy(s\left)\right)\end{array}\right. (1.1)$

例１．2 process variable derivative controller

PD_controller.slx example1_2.m

% run example1_1.m first
syms x
Px = poly2sym(num,x)/poly2sym(den,x);                               %symbolic transfer function with variable x
syms kpx kdx
Ts = Px*kpx/(1+Px*(kdx*x+kpx));                                     %close-loop transfer function with PD controller
invTs_bar = taylor(1/Ts, x, 0, 'Order', 4);                         %taylor series of first three terms
syms wn ksi
GMx = wn^2/(x^2+2*ksi*wn*x+wn^2);                                   %normal model transfer function
para_invGMs = coeffs(1/GMx, x, 'All');                              %symbolic coefficiets in matrix format
para_invTs_bar = coeffs(invTs_bar, x, 'All');                       %symbolic coefficiets in matrix format
[kpsol,kdsol] = solve(para_invGMs == para_invTs_bar, [kpx,kdx]);    %solve
kp = double(subs(kpsol, [wn ksi], [0.5 0.7]));                      %numeric substitution sym2mat
kd = double(subs(kdsol, [wn ksi], [0.5 0.7]));                      %numeric substitution sym2mat

example1_2.m shows how to calculate portion parameter $k_p$ and differential parameter $k_d$ according to specified damping rate $\zeta$ and natural frequency ${\omega }_n$ for a high order system in Fig. 1.1

Q: dynamic simulation result changing with parameters

---

→ １．2現代制御理論における高次システムの取り扱い

例１．3 state feedback control

This example shows the relevance between state feedback and process variable derivative.
The process variable derivative control in Fig. 1.1 can be rewritten as

$u\left(t\right)=\left[\begin{array}{cccc}0& 0& -k_p& -k_d\end{array}\right]{\left[\begin{array}{cccc}{x}_1\left(t\right)& x_2\left(t\right)& x_3\left(t\right)& x_4\left(t\right)\end{array}\right]}^T+k_p{y}^{ref}\left(t\right) (1.2)$

where $x\left(t\right)={\left[\begin{array}{cccc}{x}_1\left(t\right)& x_2\left(t\right)& x_3\left(t\right)& x_4\left(t\right)\end{array}\right]}^T={\left[\begin{array}{cccc}{z}_1\left(t\right)& {\dot{z}}_1\left(t\right)& y\left(t\right)& \dot{y}\left(t\right)\end{array}\right]}^T$ in Eq. 1.1.
When taking account of $z_1\left(t\right)$ , a new process variable derivative controller may be proposed as

$\begin{array}{c}u\left(t\right)=k_{p1}\left(y^{ref}\right(t)-z_1(t\left)\right)-k_{d1}{\dot{z}}_1\left(t\right)+k_{p2}\left(y^{ref}\right(t)-y(t\left)\right)-k_{d1}\dot{y}\left(t\right)\\ \iff u\left(t\right)=\left[\begin{array}{cccc}-k_{p1}& -k_{d1}& -k_{p2}& -k_{d2}\end{array}\right]x\left(t\right)+(k_{p1}+k_{p2})y^{ref}\left(t\right)\end{array} (1.3)$

example1_3.slx example1_3.m is a straight way to calculate corresponding parameters in Fig. 1.3.

% run example1_1.m first
%% transfer function
syms kp1 kp2 kd1 kd2 yref
us2 = kp1*(yref-cancelZ)-kd1*s*cancelZ+kp2*(yref-y)-kd2*s*y;        %control input
close_loop2 = us2*Ps==y;                                            %close loop
Ts2 = subs(solve(close_loop2, y),yref,1);                           %transfer function
[symnum2,symden2] = numden(Ts2);                                    %extract num den
num2 = coeffs(symnum2, s, 'All');                                   
den2 = coeffs(symden2, s, 'All');                                   %coefficients of s

ss_Ps = ss(tf);                                                     %derive A,B,C matrix for simulation
                        %WARNING:ss_Ps is not suitable for
                        %simulation 1_3 since states change,
                                                                 
k = [-1.9, -1.71, -5.84, -4.45];                                    %controller gains from book
h = 7.75;

A = [0 1 0 0;-4 -2 4 2;0 0 0 1;2 1 -2 -1];                          %A,B,C matrix by input
B = [0;2;0;0];
C = [0,0,1,0];

E: ss_Ps in example1_3.m mustn't be used in example1_3.slx as State-Space shown in Fig. 1.3 , since states change in the process of tf2ss.
In order to set state-space parameters correctly, example1_3_sup.m is proposed.

% this supplyment shows how to generate state space from differential equations
syms x1(t) x2(t) x3(t) x4(t) u;

eqs = [diff(x1(t),t) == x2(t),...
       0.5*diff(x2(t),t) == -2*(x1(t)-x3(t))-(x2(t)-x4(t))+u,...
       diff(x3(t),t) == x4(t),...
       diff(x4(t),t) == 2*(x1(t)-x3(t))+(x2(t)-x4(t))];
[vf,Yvar] = odeToVectorField(eqs);                                  % ODE to vector field

rearr = matlabFunction(Yvar);
x_state_rearr = rearr(x1,x2,x3,x4);
x_state = [x1;x2;x3;x4];
[i,j] = find(x_state*transpose(1./x_state_rearr)==1);
TrM = zeros(4);
TrM(sub2ind([4,4],i,j)) = 1;                                        %since odeToVectorField rearrange states, this
                        %step gets transformation matrix to arrange 
                        %correctly
sort_vf = TrM*vf;
F_vf = matlabFunction(sort_vf);
x_vf = F_vf(Yvar,u);
mat_A = double(equationsToMatrix(x_vf ==0, TrM*Yvar));
mat_B = double(equationsToMatrix(x_vf ==0, u));

mat_C = [0 0 1 0];                                                  %use new syms and equationsToMatrix when complicated

mat_A , mat_B and mat_C are the corresponding parameters in Fig. 1.3 . The controller parameters $k=\left[\begin{array}{cccc}-k_{p1}& -k_{d1}& -k_{p2}& -k_{d2}\end{array}\right]$ and $h=k_{p1}+k_{p2}$ in Fig. 1.3 and Eq. 1.3 may be derived by minimizing cost function

$J={\int }_0^{\infty }\left\{60\left(y\right(t)-y_c^{ref}{)}^2+u(t{)}^2\right\}dt (1.4)$

According to N: optimal regulator theory , the result of optimization Eq. 1.4 is

$k=\left[\begin{array}{cccc}-1.9& -1.71& -5.84& -4.45\end{array}\right], h=7.75 (1.5)$

Q: How to derive the result of Eq. 1.5 .
The simulation results of Fig. 1.3 with Eq. 1.5 and mat_A , mat_B , mat_C in example1_3_sup.m are shown below. A control group is established based on state-space form of Fig. 1.1
example1_3_control.slx where the corresponding state space parameters are the same with Fig. 1.4 given by example1_3_sup.m , controller gain $k_d, k_p$ are worked out by example1_2.m . The simulation result of states are shown below. Conclusively, comparing Fig. 1.2 to Fig. 1.4a and Fig. 1.6 to Fig. 1.4b , the later state-space method, which uses optimal regulator theory, is much more advantageous than previous traditional method, which exploits damping rate and natural frequency of transfer function's taylor series.
We also present another control group of common PD controller. The brief diagram is shown below. normal_pd.m

% run example1_1.m first
syms kpp kdd yref
eq3 = y == Ps*kpp*(yref-y)+kdd*s*(yref-y);
tf_pd = subs(solve(eq3, y),yref,1);
[sym_pd_num,sym_pd_den] = numden(tf_pd);
pd_num = coeffs(sym_pd_num, s, 'All');
pd_den = coeffs(sym_pd_den, s, 'All');

tl4 = taylor(1/tf_pd, s, 0, 'Order', 4);
tl5 = taylor(1/tf_pd, s, 0, 'Order', 5);
tl6 = taylor(1/tf_pd, s, 0, 'Order', 6);

By calculating 3rd to 6th order taylor series as tl4, tl5 and tl6 in normal_pd.m , the result shows that derivative gain $k_d$ only appears in taylor series with order higher than 5, which means derivative gain in common PD controller has little influence on such high order target plant with transfer function Ps in example1_1.m .

---

→ １．3多入力多出力システムに対する古典制御理論の限界

例１.4 2 inputs 2 outputs transfer function representation

example1_4.m

%% MIMO state space
syms x1(t) x2(t) x3(t) x4(t) u1 u2;

eqs = [diff(x1(t),t) == x2(t),...
       0.5*diff(x2(t),t) == -2*(x1(t)-x3(t))-(x2(t)-x4(t))+u1,...
       diff(x3(t),t) == x4(t),...
       diff(x4(t),t) == 2*(x1(t)-x3(t))+(x2(t)-x4(t))+u2];
[vf,Yvar] = odeToVectorField(eqs);
rearr = matlabFunction(Yvar);
x_state_rearr = rearr(x1,x2,x3,x4);
x_state = [x1;x2;x3;x4];
[i,j] = find(x_state*transpose(1./x_state_rearr)==1);
TrM = zeros(4);
TrM(sub2ind([4,4],i,j)) = 1;

sort_vf = TrM*vf;
F_vf = matlabFunction(sort_vf);
x_vf = F_vf(Yvar,u1,u2);
mat_A = double(equationsToMatrix(x_vf ==0, TrM*Yvar));
mat_B = double(equationsToMatrix(x_vf ==0, [u1,u2]));
mat_C = [1 0 0 0; 0 0 1 0];

%% MIMO transfer function
syms y1 y2 s;
eq1 = 0.5*s^2*y1==u1-2*(y1-y2)-(s*y1-s*y2);
eq2 = s^2*y2==u2+2*(y1-y2)+(s*y1-s*y2);
cancely1 = solve(eq2, y1);
uy1 = subs(eq1, y1, cancely1);
tf_y1u1 = subs(solve(uy1, y2),[u1,u2],[1,0]);
tf_y1u2 = subs(solve(uy1, y2),[u1,u2],[0,1]);

cancely2 = solve(eq2, y2);
uy2 = subs(eq1, y2, cancely2);
tf_y2u1 = subs(solve(uy2, y1),[u1,u2],[1,0]);
tf_y2u2 = subs(solve(uy2, y1),[u1,u2],[0,1]);

mat_tf = [tf_y1u1, tf_y1u2; tf_y2u1, tf_y2u2];

Control of an MIMO system based on transfer function representation requires an approximation to decouple the transfer function matrix, even though, there still exists interference between inputs and outputs. While control method based on state-space representation has no such worry.
Q: optimal regulator theory solved example 1.6