latest lecture

### [lecture #1] 2018.9.6 outline of the lecture, review of classical and modern control theory (1/3) †

• review : stabilization of SISO unstable plant by classical and modern control theory
• transfer functions / differential equations
• poles / eigenvalues
• impulse response / initial value response
• ...
%-- 2018/09/06 12:48 --%
s = tf('s')
P = 1/(s-1)
pole(P)
impulse(P)
k = 2
Tyr = feedback(P*k, 1)
help feedback
Tyr = feedback(P*k, 1)
step(Tyr)
k
k = 10
Tyr = feedback(P*k, 1)
step(Tyr)
• Q: Q3とQ4で話した内容で
古典制御：初期値が0である代わりにインパルスを入れる
現代制御：初期値がある代わりにインパルスを入れない
であっていますか？正確に聞きそびれたので教えてください
• A: 合っています。H∞制御の導入のところでまた詳しく説明する予定です。

### [lecture #2] 2018.9.20 review of classical and modern control theory (2/3) with introduction of Matlab/Simulink †

1. introduction of Matlab and Simulink text_fixed.pdf Basic usage of MATLAB and Simulink used for 情報処理演習及び考究II/Consideration and Practice of Information Processing II: Advanced Course of MATLAB
• interactive system (no compilation, no variable difinition)
• m file
1. system representation: Transfer Function(TF) / State-Space Representation (SSR)
• example: mass-spring-damper system
• difinition of SSR
• from SSR to TF
• from TF to SSR: controllable canonical form
2. open-loop characteristic
• open-loop stability: poles and eigenvalues
• Bode plot and frequency response ex0920_1.m mod0920_1.mdl
• cut off frequency; DC gain; -40dB/dec; variation of c
• relation between P(jw) and steady-state response
3. closed-loop stability
• Nyquist stability criterion (for L(s):stable)
• Nyquist plot ex0920_2.m mod0920_2.mdl
• Gain Margin(GM); Phase Margin(PM)
%-- 2018/09/20 13:05 --%
demo
t = [1, 2, 3]
who
a = t'
who
a(5) = 1
a(5,4) = 100
pwd
ex0920_1
who
P
P.den
P.den{:}
P
eig(P)
pole(P)
Pss = ss(P)
Pss
Pss.a
eig(Pss.a)
ex0920_2

### [lecture #3] 2018.9.27 review of classical and modern control theory (3/3) †

1. LQR problem
• controllability
• cost function J >= 0
• (semi)-positive definiteness
2. solution of LQR problem
• closed loop stability ... Lyapunov criterion
• Jmin lqr.pdfproof4.pdf (from B3「動的システムの解析と制御」)
3. example ex0927_1.m mod0927_1.mdl
%-- 2018/09/27 13:01 --%
ex0927_1
A
B
Uc
det(Uc)
eig(Q)
P
A'*P
A'*P + P*A + Q - P*B/R*B'*P
A'*P + P*A + Q - P*B*inv(R)*B'*P
Jmin
J
plot(t, J)
• Q: 計算する時間もっとほしいです
• A: もう少し見回って進捗を見るようにします。すみません。

### [lecture #4] 2018.10.4 relation between LQR and H infinity control problem (1/2) †

• GOAL: to learn difference in concepts between LQR problem and H infinity control problem
1. a simple example relating LQR and H infinity control problems
• For given plant G $G = \left[\begin{array}{c|c:c} a & 1 & b \\ \hline \sqrt{q} & 0 & 0 \\ 0 & 0 & \sqrt{r} \\ \hdashline 1 & 0 & 0 \end{array} \right] = \left\{ \begin{array}{l} \dot x = ax + bu + w\\ z = \left[ \begin{array}{c} \sqrt{q} x \\ \sqrt{r} u \end{array}\right] \\ x = x \end{array}\right.$ with zero initial state value x(0) = 0, find a state-feedback controller $u = -f x$ such that \begin{eqnarray} (i) &&\quad \mbox{closed loop is stable} \\ (ii) &&\quad \mbox{minimize} \left\{\begin{array}{l} \| z \|_2 \mbox{ for } w(t) = \delta(t) \quad \mbox{(LQR)} \\ \| T_{zw} \|_\infty \mbox{($H_\infty$ control problem)}\end{array}\right. \end{eqnarray}
• comparison of norms in (ii) (for a = -1, b = 1, q = 1, r = 1) $\begin{array}{|c||c|c|}\hline & \mbox{LQR}: f=-1+\sqrt{2} & \quad \quad H_\infty: f=1\quad\quad \\ \hline\hline J=\|z\|_2^2 & & \\ \hline \|T_{zw}\|_\infty & & \\ \hline \end{array}$
2. an alternative description to LQR problem
1. J = (L2 norm of z)^2
2. impulse resp. with zero initial value = initial value resp. with zero disturbance
3. definition of H infinity norm (SISO)
s = tf('s');
G1 = 1/(s+1);
bode(G1);
norm(G1, 'inf')
G2 = 1/(s^2 + 0.1*s + 1);
bode(G2);
norm(G2, 'inf')
4. definition of H infinity norm (SIMO)
5. solve the problem by hand
6. solve the problem by tool(hinfsyn) ex1004.m
%-- 2018/10/04 12:58 --%
s = tf('s');
G1 = 1/(s+1);
bode(G1);
norm(G1, 'inf')
G2 = 1/(s^2 + 0.1*s + 1);
bode(G2);
norm(G2, 'inf')
format long e
norm(G2, 'inf')
grid on
ex1004
• Q: z を L して Z(s) = Tzw W(s) の Tzw の導出？ についていけなかった。
• A: その直前の話 $X(s) = \frac{1}{s - (a - bf)}W(s) = \frac{1}{s + 1 + f}W(s)$ が大丈夫とすると、 $z = \left[ \begin{array}{c} \sqrt{q} \\ -\sqrt{r} f \end{array}\right] x= \left[ \begin{array}{c} 1 \\ -f \end{array}\right]x$ より、 $Z(s) = \left[ \begin{array}{c} 1 \\ -f \end{array}\right] \frac{1}{s + 1 + f}W(s)$ すなわち、 $T_{zw}(s) = \left[ \begin{array}{c} 1 \\ -f \end{array}\right] \frac{1}{s + 1 + f} = \left[ \begin{array}{c} \frac{1}{s + 1 + f} \\ \frac{-f}{s + 1 + f} \end{array} \right]$ です。

### [lecture #5] 2018.10.11 relation between LQR and H infinity control problem (2/2) †

1. complete the table in simple example
2. confirm the cost function J for both controllers by simulation mod1011.mdl
• block diagram in the simulink model
• how to approximate impulse disturbance
• impulse disturbance resp. with zero initial condition = initial condition resp. with zero disturbance
3. confirm the closed-loop H infinity norm for both controllers by simulation
• H infinity norm = L2 induced norm
• review: steady-state response; the worst-case disturbance w(t) which maximizes L2 norm of z(t) ?
• what is the worst-case disturbance in the simple example ?
4. general state-feedback case: hinf.pdf
• includes the simple example as a special case
• LQR lqr.pdf is included as a special case in which gamma -> infinity, w(t) = 0, B2 -> B, and non-zero x(0) are considered
%-- 2018/10/11 12:58 --%
sqrt(2-sqrt(2))
mod1011
f
ex1004
a
b
q
r
A
f
f = -1+sqrt(2)
x0
x0 = 0
h
h = 0.1
zz
h = 0.001
zz
f = 1
zz
x0 = 1
zz
f
h
h = 1
x0
x0 = 0
zz
zz(end)
zz(end)/ww(end)
sqrt(zz(end)/ww(end))
h
h = 100
sqrt(zz(end)/ww(end))
f = sqrt(2)-1
h
sqrt(zz(end)/ww(end))
• Q: LQR は外乱を考慮しない制御理論という認識でよいか
• A: 元の問題に外乱の項が無い、という意味ではそうです。が、初期値応答を考慮することとインパルス外乱を考慮することは等価なので、その意味では外乱を考慮していると言えます。 一方、明示的に外乱（ノイズ）の項を加えたLQG問題もあります。

### [lecture #6] 2018.10.18 Mixed sensitivity problem 1/3 †

1. outline: map_v1.1_mixedsens1.pdf
• sensitivity function S and complementary sensitivity function T
2. H infinity control problem (general case)
• with generalized plant G
• including the state-feedback case)
3. reference tracking problem
• corresponding generalized plant G ?
• how to translate the condition (ii) into one with H infinity norm ?
• introduction of weighting function for sensitivity function in (ii)
4. design example ex1018_1.m ex1018_2.m
5. the small gain theorem
• proof: Nyquist stability criterion
%-- 2018/10/18 12:58 --%
ex1018_1
P
pole(P)
eig(P)
ex1018_2
G
who
K_hinf
eig(K.a)
eig(K_hinf.a)
ex1018_2
K_hinf
figure(5)
bode(K_hinf)
• Q: ${}^\forall \ \omega$ とか ${}^\forall \ \Delta$ とかの ∀って何でしたっけ。
• A: 「すべての」です。
• Q: ex1018_2.m の 13行目の説明がはやくて分からなかった。
• A: hinfsyn ... H∞制御問題を解きます。第一引数に一般化プラント（今の例ではG）を指定します。第二、第三引数に、制御入力u と 観測出力yのサイズ（今の例ではスカラ）を指定します。 戻り値については、一つ目から順にH∞ノルムγを最小化するコントローラ（K_hinf）、そのときの閉ループ系（clp）、γ（gopt）、です。その場ですぐ聞いてもらえると他の人も助かると思います。

### [lecture #7] 2018.10.25 Mixed sensitivity problem 2/3 †

1. review: the small gain theorem ... robust stability = H infinity norm condition
2. outline: from point to set map_v1.1_mixedsens2.pdf
3. normalized uncertainty Delta
4. uncertainty model
5. how to determine P0 and WT
• example: frequency response of plant with perturbation ex1025_1.m
• frequency response based procedure for P0 and WT ex1025_2.m
6. robust stabilization problem and equivalent problem
%-- 2018/10/25 13:25 --%
ex1025_1
ex1025_2
ex1025_3
mod1025
c
c = 0.8
c = 1.3
c = 2
• Q: additive uncertainty model は multiplicative 〜 を含むということですか？
• A: 確かに今日説明したように、W(s) = P0(s)WT(s) とおけば、W(s)は加法的摂動モデルの重み関数となります。その意味では形式的に、含みます。ただその場合は、W(s) がノミナルプラントと同じ動特性を持つことになりW(s)の次数がムダに高くなるため、加法的摂動モデルではなく乗法的摂動モデルをそのままその後の設計問題に用います。つまり、実際にコントローラの設計を行う上では、加法的摂動モデルは乗法的摂動モデルを含みません。答になっていない場合はまた聞いてください。
• Q:
• Q: ex1025_1.m で $\tilde P(s) = \frac{1}{(s-1)(cs+1)}$ の (cs+1) の部分が不安定にならない理由がよく分かりませんでした
• A: その部分に対応する極は cs + 1 = 0 より、s = -1/c です。今、c は正の実数として与えられているので、-1/c は負、つまり安定です。

### [lecture #8] 2018.11.1 Mixed sensitivity problem 3/3 †

• review: map_v1.1_mixedsens2.pdf (1)robust stabilization and (2)performance optimization
• outline:
1. how to design controllers considering both conditions in (1) and (2)
2. gap between nominal performance and robust performance
1. mixed sensitivity problem ---> (1) and (2) : proof
2. generalized plant for mixed senstivity problem
3. design example ex1101_1.m minimize gamma by hand
4. gamma iteration by bisection method ex1101_2.m
5. nominal performance and robust performance ex1101_3.m
6. introduction of robust performance problem
%-- 2018/11/01 12:59 --%
ex1101_1
gam
ex1101_1
K
ex1101_2
ex1101_3

### [lecture #9] 2018.11.8 robust performance problem 1/3 †

1. review
• mixed sensitivity problem : N.P. but not R.P.
• robust performance problem (R.P.) c.f. the last whiteboard
2. an equivalent robust stability (R.S.) problem to R.P.
• with structured uncertainty Delta hat
3. definition of H infinity norm for general case (MIMO)
• definition of singular values and the maximum singular value
M = [1i, 1i; 0, 1]
M'
eig(M'*M)
svd(M)
• mini report #1 report1.pdf ... You will have a mini exam #1 related to this report
4. proof of ||Delta hat||_inf <= 1
5. design example: ex1108_1.m
• robust performance is achieved but large gap
• non structured uncertainty is considered ... the design problem is too conservative
%-- 2018/11/08 13:03 --%
ex1101_1
ex1101_2
ex1101_3
1i
j
M = [1i, 1i; 0, 1]
M'
eig(M'*M)
svd(M)
(3+sqrt(5))/2
sqrt((3+sqrt(5))/2)
ex1108_1
• Q: 共役複素数の~（バー）と最大を表す~（バー）は行列がR（実数）かC（複素数）かで使い分けされているということですか。
• A: そういうことではなく、単に、複素共役にも最大にも同じ記号を当てています。確かに言われてみると紛らわしいですが、文脈から混同する恐れはないので、このままにさせてください。
• Q: λ(M*M)が実数の証明レポートにしません？
• A: 「mini」レポートにならなさそうなので、今年度はレポートにはしません。良いアイデアをありがとうございます。

### [lecture #10] 2018.11.15 Robust performance problem (2/3) †

1. return of mini report #1
2. review
• robust performance but too conservative
ex1108_1
• robust stability problem for Delta hat and its equivalent problem(?) with Delta tilde
• structured unertainty Delta hat and unstructured uncertainty Delta tilde
3. SVD: singular value decomposition
• definition
• meaning of the largest singular value
• 2-norm of vectors (Euclidean norm)
• SVD for 2-by-2 real matrix ex1115_1.m
%-- 2018/11/15 13:01 --%
M = [1, 2; 3, 4]
[U, S, V] = svd(M)
M = [1, 1i/sqrt(2); 1, -1i/sqrt(2)]
[U, S, V] = svd(M)
U'*U
format long e
U'*U
V'*V
M = [1, 2; 3, 4]
[U, S, V] = svd(M)
fomat short
format short
[U, S, V] = svd(M)
ex1115_1
S
sys = tf(1, [1, 1])
sys2 = tf(1, [1, 1])
sys2 = tf(1, [1, 2])
sys = tf(1, [1, 1])
sys
sys2
M = [0, sys; 0, sys2]
M
hinfnorm(M)
M = [0, sys; 0, sys]
hinfnorm(M)
• Q: $\sigma_i u_i = M v_i$ が証明できれば $\bar \sigma(M) = \max_{a\neq 0}\frac{\| b\|_2}{\|a\|_2}$ がいえるということで合っていますか？
• A: 合っています。（厳密には、ベクトル a を $a = \Sigma_{i=1}^m \alpha_i v_i$ のように、ベクトル vi の線形結合で表して実数 αi を振って2ノルムの比を最大化すると、結局 a = v1 （ただし v1 は最大特異値に対応する V のベクトル）となります。
• Q: $\|M\|_\infty < 1$ をたしかめる時 $\bar \sigma(M)$ を計算して出しますが 1より大きいかどうか考えるなら最後の√をつけなくても良くないですか？
• A: その通りです。
• Q: Σは M の大きさによっては、 $\Sigma = \left[ \begin{array}{c} \begin{array}{ccc} \sigma_1 & & 0 \\ & \ddots & \\ 0 & & \sigma_m \end{array} \\ \hdashline 0 \mbox{←ここ} \end{array} \right]$ のゼロ行列はなくなりますか？ また、 $M^{p\times m}$ が p < m のときは Σ はどうなりますか？
• A: Mが正方行列つまりp=mのときゼロ行列はなくなります。また、p<mのときはΣは横長になり、0行列は下でなくて右につきます。

### [lecture #11] 2018.11.22 Robust performance problem (3/3) †

1. review: R.S. prob. for Delta hat and Delta tilde
2. scaled H infinity control problem
3. how to determine structure of scaling matrix
4. design example ex1122_1.m
ex1108_1
gam2 = gam_opt
K2 = K_opt;
ex1122_1
gam_opt
5. mini report #2 report2.pdf
6. mini exam #1 (10 min.)
%-- 2018/11/22 13:37 --%
ex1108_1
gam2 = gam_opt
K2 = K_opt
ex1122_1
gam_opt
gam2
d_opt

### [lecture #12] 2018.11.29 Robust performance problem (3/3) (cont.), Control system design for practical system (1/3) †

1. return of mini exam #1 and mini report #2
2. review of scaling ex1129_1.m
3. introduction of a practical system: Speed control of two inertia system with servo motor
• experimental setup
setup.pdf ... spring plate'' in Table 1 is replaced with a thin long shaft
photo1.png
• objective of control system = disturbance attenuation control problem: to drive the drive-side servomotor by generating proper driving signal u using drive-side speed (as well as driven-side speed) y such that the driven-side speed is maintained at constant against the disturbance torque input w
• frequency response experiment and physical model of speed control system
• room 374 @ Dept. Mech. Bldg. 2
%-- 2018/11/29 13:07 --%
ex1129_1
gam2
gam3
clp2.b
format long
norm(M3_d, 'inf')
ex1129_2

### [lecture #13] 2018.12.6 Control system design for practical system (2/3) †

1. mini exam #2
2. review of the experimental system
• closed-loop system of 2-by-2 plant G and controller K
• closed-loop gain is desired to be minimized
• how to handle modeling error of G ?
3. design example
• frequency response experiment data
servo1.dat
servo2.dat
• determination of plant model(nominal plant and additive uncertainty weight)
nominal.m
weight.m
• configuration of generalized plant and controller design by scaled H infinity control problem using one-dimensional search on the scaling d
cont.m
• comparison of closed-loop gain characteristics with and without control
compare.m
• result of control experiment
result.dat
perf.m
4. final report and remote experimental system
1. design your controller(s) so that the system performance is improved compared with the design example
2. Draw the following figures and explain the difference between two control systems (your controller and the design example):
1. bode diagram of controllers
2. gain characteristic of closed-loop system from w to z
3. time response of control experiment
3. Why is the performance of your system improved(or unfortunately deteriorated)?
• due date: 4th(Fri) Jan 17:00
• submit your report(pdf or doc) by e-mail to kobayasi@nagaokaut.ac.jp
• You can use Japanese
• maximum controller order is 20
• the system will be started until next lecture
• You can send up to 510 controllers
• control experimental results will be uploaded here
• freqresp ... frequency response will be measured and uploaded everyday
5. how to improve the performance ?
• order of the nominal plant accuracy of the nominal(physical) model
• weighting for robust stability
6. specifications of the experimental system
1. program sources for frequency response experiment
• freqresp.h
• freqresp_module.c
• freqresp_app.c
• format of servo1.dat (w is used instead of u for servo2.dat)
1st column ... frequency (Hz)
2nd column ... gain from u(Nm) to y(rad/s)
3rd column ... phase (deg) from u to y
4th column ... gain from u to z
5th column ... phase (deg) from u to z
2. program sources for control experiment
• hinf.h
• hinf_module.c
• hinf_app.c
• format of result.dat
1st column: time (s)
4th column: u (Nm)
5th column: w (Nm)
3. configuration of control experiment
• disturbance signal w is specified as described in hinf.h and hinf_module.c:
w = 0; // disturbance torque for driven motor
if((t > 2)&&(t < 3)){
w = RATED_TORQ * -0.15;
}
if((t > 4)&&(t < 5)){
w = RATED_TORQ * -0.1 * sin(2*M_PI*5.0 * (t-4.0));
}
da_conv(torq_volt_conv_1(w), 1);
• control signal u is limited as specified in hinf.h and hinf_module.c:
#define U_MAX (RATED_TORQ / 3.0)

if(u > U_MAX) u = U_MAX;
if(u < -U_MAX) u = -U_MAX;
u is generated by PI control for t < 1(s). Your designed controller is started at t = 1(s).
4. calculation of rotational speed
• The rotational speed is approximately calculated by using difference for one sampling period in hinf_module.c and freqresp_module.c like:
theta_rad[0] = (double)read_theta(0) / (double)Pn212 * 2.0 * M_PI;
theta_rad_before[1] = theta_rad[1];
where the sampling period is given as 0.25 ms.
%-- 2018/12/06 12:53 --%
nominal
w
size(w)
servo1
size(servo1)
nupwd
pwd
cd data
ls
cd ..
nominal
size(w)
size(resp)
size(servo2)
size(servo1)
n
nominal
n
k
resp
k
w
clear all
nominal
weight
cont
compare
perf
nominal
• Q: コントローラは何で書くのか？
• A: 説明が曖昧でしたが、コントローラの状態空間実現（のファイル）を作るだけで良いです。 具体的には、これらのファイルが実験機で動作しているC言語のプログラムにロードされ、制御実験が行われます。

### [lecture #14] 2018.12.13 Control system design for practical system (3/3) †

• return of mini exam #2
• review & supplemental explanations
• final report
• control objective is to suppress speed fluctuation of the driven-side motor not the drive-side motor (a large vibration at the drive-side motor is allowed)
• no strict control objective is given ( there is a freedom to define what is good performance; a frequency dependent weighting function can be introduced to evaluate the performance )
• c2d() is used to discretize the resultant continuous-time controller
• web based remote experiment system
• now you can login after registration
• room temperature is displayed and stored in temp.txt (Bosch Sensortec BME280 is used)
• the system will be unavailable from 21(Fri) 7:30 on 22(Sat) to 19:00 on 22(Sat) due to electrical construction scheduled on 22(sat)
• preparation of your own controller(s) by using the remote experiment system
%-- 2018/12/13 13:58 --%
help subspace
nominal_ave
pwd
cd ..
pwd
nominal_ave
G_g
G1_g
G2_g
(G_g + G10_g)/2
bode((G_g + G10_g)/2)
figure(2)
bode((G_g + G10_g)/2)
figure(2)

### [lecture #15] 2018.12.20 Control system design for practical system (cont.) †

• preparation of your own controller(s)
• use design example #3 for comparison with your design
• discuss relationship between the required figures (a), (b), and (c)
• driving and driven torque are set to zero when driven-motor speed z exceeds a limit as in hinf.h:
#define SPEED_MAX 60.0 // driving torque becomes zero when rotational speed is out of range from -60.0 rad/s to 60.0 rad/s in order to avoid shutting down the driving motor by the maximum speed excess alarm (2018.12.19)
and hinf_module.c:
       int flag_speed_excess = 0; // 2018.12.19

if(flag_speed_excess == 1) w = 0; // 2018.12.20
da_conv(torq_volt_conv_1(w), 1);

if(fabs(z) > SPEED_MAX) flag_speed_excess = 1; // 2018.12.20

if(flag_speed_excess == 1) u = 0; // 2018.12.19

#ifndef NO_CONTROL
da_conv(torq_volt_conv_0(u), 0);
#endif
• the system will be unavailable 7:30 -- 19:00 on 22(Sat)
• questionnaires
• to university
• for web-based experimental environment
%-- 2018/12/20 13:19 --%
nominal
weight
cont
compare
nominal
weight
ctrlpref
cont
compare
weigh
weight
cont
weight
nominal
weight
• Q: うまくいくとやっぱりたのしいですね。
• A: モチベーションが高まったところで、過去の授業の内容を復習すると理解が深まりなお良いと思います。
• 2018.12.20 19:00 due to unexpected alert occurred at 17:33, hinf_module.c is further modified at 18:40 as follows (drive-side speed is additionally evaluated to prevent the alert):
	    //	    if(fabs(z) > SPEED_MAX) flag_speed_excess = 1; // 2018.12.20
if((fabs(z) > SPEED_MAX) || (fabs(y) > SPEED_MAX)) flag_speed_excess = 1; // 2018.12.20 further modified after today's lecture

Last-modified: 2018-12-20 (木) 19:03:24 (247d)