授業
Advanced Automation†
[lecture #1] 2016.9.1 outline of the lecture, review of classical and modern control theory (1/3)†
- outline of this lecture
- syllabus
- evaluation
- mini report #1 ... 10%
- mini exam #1 ... 10%
- mini report #2 ... 10%
- mini exam #2 ... 10%
- final report ... 60%
- schedule2016 (tentative)
- map
- review : stabilization of SISO 1st-order unstable plant by classical and modern control theory
- transfer functions / differential equations
- poles / eigenvalues
- impulse response / initial condition response
- ...
%-- 2016/09/01 14:26 --%
s = tf('s')
k = 2
alpha = -1
Tyr = k/(s+alpha+k)
step(Tyr, 'b')
k = 100
Tyr2 = k/(s+alpha+k)
step(Tyr, 'b', Try2, 'r--')
step(Tyr, 'b', Tyr2, 'r--')
[lecture #2] 2016.9.8 review of classical and modern control theory (2/3) with introduction of Matlab/Simulink†
- 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
- 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
- open-loop characteristic
- open-loop stability: poles and eigenvalues
- Bode plot and frequency response ex0908_1.m mod0908_1.mdl
- cut off frequency; DC gain; -40dB/dec; variation of c
- relation between P(jw) and steady-state response
- closed-loop stability
- Nyquist stability criterion (for L(s):stable)
- Nyquist plot ex0908_2.m mod0908_2.mdl
- Gain Margin(GM); Phase Margin(PM)
%-- 2016/09/08 13:08 --%
a = 1
a
pwd
ls
ex0908_1
sqrt(m/k)/(2*pi)
sqrt(k/m)/(2*pi)
[lecture #3] 2016.9.15 review of classical and modern control theory (3/3)†
- LQR problem
- controllability
- cost function J >= 0
- (semi)-positive definiteness
- solution of LQR problem
- ARE and quadratic equation
- closed loop stability ... Lyapunov criterion
- Jmin
lqr.pdf ≒ proof4.pdf (from B3「動的システムの解析と制御」)
- example
mod0915.mdl
A = [1, 2; 0, -1]; % unstable plant
B = [0; 1];
Uc = ctrb(A,B);
det(Uc) % should be nonzero
C = eye(2); % dummy
D = zeros(2,1); % dummy
F = [0, 0]; % without control
x0 = [1; 1]; % initial state
Q = eye(2);
R = 1;
P = are(A, B/R*B', Q);
eig(P) % should be positive
F = R\B'*P;
x0'*P*x0
%-- 2016/09/15 14:18 --%
mod0915
A = [1, 2; 0, -1]; % unstable plant
B = [0; 1];
Uc = ctrb(A,B);
det(Uc) % should be nonzero
C = eye(2); % dummy
D = zeros(2,1); % dummy
F = [0, 0]; % without control
x0 = [1; 1]; % initial state
Q = eye(2);
R = 1;
P = are(A, B/R*B', Q);
eig(P) % should be positive
F = R\B'*P;
x0'*P*x0
x0
J
x0'*P*x0
plot(t, J)
- Q: 各行列の意味がよくわからなかった。
- A: Matlab で入力した行列のことと仮定して回答します(そうでない場合はまた聞いてください)。
A = [1, 2; 0, -1]; % unstable plant
B = [0; 1];
は、制御対象の例として与えました。何でも良いですが、ここでは不安定な方が分かり易いかと思い、不安定な制御対象としました。
Uc = ctrb(A,B);
det(Uc) % should be nonzero
は、可制御性を確認している部分です。(A,B)の可制御性行列を Uc として作り、det で行列式を表示しています。
C = eye(2); % dummy
D = zeros(2,1); % dummy
おそらくこの部分を全く説明しなかったので混乱していると思います(すみません)が、
この C と D を y = Cx + Du に代入すると、y = x となります。つまり、出力に状態をそのまま取り出すための設定です。これで、Simulink の制御対象のブロックの出力が、状態ベクトル x となります。
F = [0, 0]; % without control
フィードバック制御無しの設定です。ここまで入力して一度、シミュレーションスタートするべきでした(不安定な、信号が発散する様子が確認できます)。が、時間なくスキップしたと思います。すみません。
x0 = [1; 1]; % initial state
初期値ベクトルを適当に与えています。
Q = eye(2);
R = 1;
評価関数の重み行列(という呼び方で説明しなかったと思いますが、そう呼ばれます)の設定です。制御系の設計者が与えるものです。
P = are(A, B/R*B', Q);
リカッチ方程式の正定解Pを計算します。この結果を実際にリカッチ方程式に代入すると0となるところも見せるべきでした。やってみてください。
eig(P) % should be positive
全て正の実数となることから、実際に正定となっていることを確認します。
F = R\B'*P;
状態フィードバックゲインFを計算しています。
x0'*P*x0
評価関数Jの最小値を表示します。スカラの簡単な場合について、求めた最適解fをJの式に代入しても、この同じ式になる、ということも確認できるのですが、時間なくスキップしました。やってみてください。
[lecture #4] 2016.9.29 relation between LQR and H infinity control problem (1/2)†
- GOAL: to learn difference in concepts between LQR problem and H infinity control problem
- 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 condition x(0) = 0,
find a state-feedback controller
\[ u = -f x \]
such that
\begin{eqnarray}
(i) &&\quad \mbox{closed loop is stable} \\
(ii) &&\quad
\left\{\begin{array}{l} \| z \|_2 \rightarrow \mbox{min for } w(t) = \delta(t) \quad \mbox{(LQR)} \\
\| T_{zw} \|_\infty \rightarrow \mbox{min} \quad \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}
\]
- an alternative description to LQR problem
- J and L2 norm of z
- impulse disturbance resp. with zero initial condition and initial condition resp. with zero disturbance
- definition of H infinity norm (SISO)
s = tf('s');
P1 = 1/(s+1);
bode(P1);
norm(P1, 'inf')
P2 = 1/(s^2 + 0.1*s + 1);
bode(P2);
norm(P2, 'inf')
- definition of H infinity norm (SIMO)
- solve the problem by hand
- solve the problem by tool(hinfsyn)
ex0929.m
%-- 2016/09/29 13:46 --%
s = tf('s');
P1 = 1/(s+1);
bode(P1);
norm(P1, 'inf')
P2 = 1/(s^2 + 0.1*s + 1);
bode(P2);
norm(P2, 'inf')
format long e
norm(P2, 'inf')
grid on
ex0929
- Q: H∞ノルムが何を表しているのか
- A: 簡単に言うと安定なシステムの最大ゲインです。
[lecture #5] 2016.10.6 relation between LQR and H infinity control problem (2/2)†
- complete the table in simple example
- confirm the cost function J for both controllers by simulation mod1006.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
- confirm the closed-loop H infinity norm for both controllers by simulation (common mdl file is available)
- review: steady-state response (see whiteboard#4 @ lec. #2)
- H infinity norm = L2 induced norm
- how to construct the worst-case disturbance w(t) which maximizes L2 norm of z(t) ?
- what is the worst-case disturbance in the simple example ?
- general 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
%-- 2016/10/06 13:09 --%
1/sqrt(2)
sqrt(2+sqrt(2))
sqrt(2-sqrt(2))
mod1006
ex0929
A
B
C
D
h
h = 0.01
x0 = 0
f
f = -1+sqrt(2)
zz
zz(end)
h = 0.0001
zz(end)
x0
x0 = 1
zz(end)
f
f = 1
x0
x0 = 0
zz(end)
x0
h
h = 10
zz(end)/ww(end)
sqrt(zz(end)/ww(end))
f
h = 100
sqrt(zz(end)/ww(end))
f
f = -1+sqrt(2)
sqrt(zz(end)/ww(end))
whiteboard #1 ... sorry for the missing picture in which A, B1, B2, C1, and D12 were additionally written in red color.
!!! the remaining page is under construction (the contents below are from 2015) !!!
[lecture #6] 2015.10.22 Mixed sensitivity problem 1/3†
- review &ref(): File not found: "map_v1.0_intro1.pdf" at page "授業/制御工学特論2016"; and outline
- H infinity control problem (general form)
- reference tracking problem
- weighting function for sensitivity function
- design example &ref(): File not found: "ex1022_1.m" at page "授業/制御工学特論2016"; &ref(): File not found: "ex1022_2.m" at page "授業/制御工学特論2016";
- the small gain theorem
- proof: Nyquist stability criterion
- from performance optimization to robust stabilization
%-- 10/22/2015 2:06 PM --%
ex1022_1
eig(P)
ex1022_2
#ref(): File not found: "2015.10.22-1.jpg" at page "授業/制御工学特論2016"
#ref(): File not found: "2015.10.22-2.jpg" at page "授業/制御工学特論2016"
#ref(): File not found: "2015.10.22-3.jpg" at page "授業/制御工学特論2016"
[lecture #7] 2015.10.29 Mixed sensitivity problem 2/3†
- review &ref(): File not found: "map_v1.0_intro2.pdf" at page "授業/制御工学特論2016"; and outline
- an equivalent problem of robust stabilization for reference tracking problem
- uncertainty model and normalized uncertainty Delta
- robust stabilization problem and an equivalent problem
- practical example of plant with perturbation &ref(): File not found: "ex1029_1.m" at page "授業/制御工学特論2016";
- how to determine the model &ref(): File not found: "ex1029_2.m" at page "授業/制御工学特論2016";
- design example and simulation &ref(): File not found: "ex1029_3.m" at page "授業/制御工学特論2016"; &ref(): File not found: "mod1029.mdl" at page "授業/制御工学特論2016";
%-- 10/29/2015 1:52 PM --%
ex1029_1
ex1029_2
ex1029_3
mod1029
c
c = 0.8
#ref(): File not found: "2015.10.29-1.jpg" at page "授業/制御工学特論2016"
#ref(): File not found: "2015.10.29-2.jpg" at page "授業/制御工学特論2016"
#ref(): File not found: "2015.10.29-3.jpg" at page "授業/制御工学特論2016"
[lecture #8] 2015.11.5 Mixed sensitivity problem 3/3†
- review : (1)robust stabilization and (2)performance optimization
- mixed sensitivity problem : a sufficient condition for (1) and (2)
- proof by definition of H infinity norm
- construction of the generalized plant
- design example &ref(): File not found: "ex1105_1.m" at page "授業/制御工学特論2016";
- gamma iteration by bisection method &ref(): File not found: "ex1105_2.m" at page "授業/制御工学特論2016";
- a problem of the mixed sensitivity problem: nominal performance and robust performance &ref(): File not found: "ex1105_3.m" at page "授業/制御工学特論2016";
- introduction of robust performance problem
%-- 11/5/2015 1:37 PM --%
ex1105_1
ex1105_2
gam
ex1105_2
WT
ex1105_2
ex1105_3
ex1105_2
ex1105_3
#ref(): File not found: "2015.11.05-1.jpg" at page "授業/制御工学特論2016"
#ref(): File not found: "2015.11.05-2.jpg" at page "授業/制御工学特論2016"
#ref(): File not found: "2015.11.05-3.jpg" at page "授業/制御工学特論2016"
[lecture #9] 2015.11.12 robust performance problem 1/3†
- review: robust performance problem
- an equivalent robust stability problem
- definition of H infinity norm for general case (MIMO)
- definition of (maximum) singular value
M = [j, 0; -j, 1]
M'
eig(M'*M)
svd(M)
- mini report #1 report1.pdf
- write by hand
- due date and place of submission -> see schedule2015
- check if your answer is correct or not before submission by using Matlab
- You will have a mini exam #1 related to this report
- SVD: singular value decomposition
- definition
[U,S,V] = svd(M)
M = [j, 0; -j, 1; 2, 3]
- unitary matrix and 2 norm of vectors
- a property of SVD: input-output interpretation
- illustrative example: rotation matrix &ref(): File not found: "ex1112_1.m" at page "授業/制御工学特論2016";
- H infinity norm of Delta hat
%-- 11/12/2015 1:01 PM --%
M = [j, 0; -j, 1]
M'
eig(M'*M)
svd(M)
M = [j, 0; -j, 1]
M'
eig(M'*M)
svd(M)
(3+sqrt(5))/2
sqrt((3+sqrt(5))/2)
help svd
[U,S,V] = svd(M)
U'*U
M = [j, 0; -j, 1; 2, 3]
[U,S,V] = svd(M)
V'*V
ex1112_1
#ref(): File not found: "2015.11.12-1.jpg" at page "授業/制御工学特論2016"
#ref(): File not found: "2015.11.12-2.jpg" at page "授業/制御工学特論2016"
#ref(): File not found: "2015.11.12-3.jpg" at page "授業/制御工学特論2016"
#ref(): File not found: "2015.11.12-4.jpg" at page "授業/制御工学特論2016"
- Q: Σでノルムが決まる?→固有値で
- A: 与えられた行列Mを、入力をベクトル a、出力をベクトル b とする入出力システム
\[ b = M a \]
とみなすと、出力の2ノルムの最大値は、入力の2ノルムの
\[ \bar \sigma(M) \] 倍となります(それを超えるような a はない)。 質問の意図と違っている場合はまた聞いてください。
- Q: it was too fast
- A: This might be caused by my less explanation in Japanese. I will improve this in the next lecture.
[lecture #10] 2015.11.19 Robust performance problem (2/3)†
- return of mini report #1
- review and outline: robust stability problem for Delta hat and its equivalent problem(?)
- signal vector's size is not restricted in H infinity control problem and small gain theorem
- H infinity norm of Delta hat
- design example: robust performance is achieved &ref(): File not found: "ex1119_1.m" at page "授業/制御工学特論2016";
- non structured uncertainty is considered ... the design problem is too conservative
%-- 11/19/2015 1:23 PM --%
doc hinfsyn
ex1105_2
ex1105_3
gam_opt
ex1119_1
gam_opt
svd([1/sqrt(2), 0; 1/sqrt(2), 0])
#ref(): File not found: "2015.11.19-1.jpg" at page "授業/制御工学特論2016"
#ref(): File not found: "2015.11.19-2.jpg" at page "授業/制御工学特論2016"
- Q: Delta tilder is more tractability than Delta hat, Delta hat is more solvability than Delta tilder, OK?
- A: That is right. For problems
- (P1) robust stabilization against Delta hat
- (P2) robust stabilization against Delta tilde,
P1 is more solvable than P2 because of the smaller uncertain set.
P2 is more tractable than P1 because of the ignorance of the structure.
[lecture #11] 2015.11.26 Robust performance problem (3/3)†
- review
- robust performance problem with Delta hat and conservative design problem with Delta tilde
- inclusion relation between two uncertain sets
- introduction of the scaled H infinity control problem
- how to determine structure of scaling matrix
- design example moved to next lecture
% less conservative design
ex1105_2
ex1105_3
ex1119_1
gam_opt0 = gam_opt;
K_opt0 = K_opt;
#ref(): File not found: "ex1126_1.m" at page "授業/制御工学特論2016"
- effect of scaling matrix moved to next lecture
#ref(): File not found: "ex1126_2.m" at page "授業/制御工学特論2016"
- mini exam #1
#ref(): File not found: "2015.11.26-1.jpg" at page "授業/制御工学特論2016"
#ref(): File not found: "2015.11.26-2.jpg" at page "授業/制御工学特論2016"
- Q: 薄いマーカーが見辛い
- A: すみません。次回、新しいマーカーに交換してもらいます。
[lecture #12] 2015.12.3 Robust performance problem (3/3) (cont.), Control system design for practical system (1/3)†
- return of mini exam #1; schedule of mini report #2 and exam #2
- review of the scaled H infinity control problem
- comments on mu-synthesis prolem
- design example (moved from the previous lecture)
% less conservative design
ex1105_2
ex1105_3
ex1119_1
gam_opt0 = gam_opt;
K_opt0 = K_opt;
#ref(): File not found: "ex1126_1.m" at page "授業/制御工学特論2016"
- effect of scaling matrix (moved from the previous lecture)
#ref(): File not found: "ex1126_2.m" at page "授業/制御工学特論2016"
- mini report #2 report2.pdf
- write by hand
- due date and place of submission -> see schedule2015
- check if your answer is correct or not before submission by using Matlab
- You will have a mini exam #2 related to this report
- controller design for practical system: active noise control in duct
- introduction of experimental setup
#ref(): File not found: "exp_apparatus1.jpg" at page "授業/制御工学特論2016"
#ref(): File not found: "exp_apparatus2.jpg" at page "授業/制御工学特論2016"
- objective of control system: to drive control loudspeaker by generating proper driving signal u using reference microphone output y such that the error microphone's output z is attenuated against the disturbance input w
- frequency response experiment
#ref(): File not found: "ex1203_1.m" at page "授業/制御工学特論2016"
%-- 12/3/2015 1:27 PM --%
ex1105_2
ex1105_3
ex1119_1
gam_opt0 = gam_opt;
K_opt0 = K_opt;
who
gam_opt0
ex1126_1
gam_opt
d_opt
ex1126_2
ex1203_1
#ref(): File not found: "2015.12.03-1.jpg" at page "授業/制御工学特論2016"
#ref(): File not found: "2015.12.03-2.jpg" at page "授業/制御工学特論2016"
#ref(): File not found: "2015.12.03-3.jpg" at page "授業/制御工学特論2016"
- Q: Why
\[ \mbox{(i) } \hat G(0.0111) \mbox{ is stable, and} \]
\[ \mbox{(ii) } \hat G(0.0103) \mbox{ is unstable,} \]
hold ?
- A: In the first design (i), controller K_opt0 and minimized value of gamma (gam_opt0) are obtained so that the H infinity norm of the resultant closed-loop system clp0 without connecting Delta tilde is less than 1. Therefore, the closed-loop system composed of clp0 and Delta tilde is stable according to the small gain theorem.
On the other hand, in the second design (ii), controller K_opt and minimized value of gamma (gam_opt) are obtained so that the H infinity norm of the resultant closed-loop system *with scaling* (clp_d) is less than 1. Please note that there is no guarantee for the H infinity norm of the closed-loop system *without scaling* (clp_1) to be less than 1. Indeed, we confirmed that the H infinity norm of clp_1 was larger than 1 in our example, by which the closed-loop system composed of clp_1 and Delta tilde is unstable according to the small gain theorem. However, Delta hat which has diagonal structure, can be connected to clp_1 without loosing closed-loop stability. (Maybe I didn't explain the last sentence. I'm sorry for this, if this is the reason of your question.)
If this answer is not sufficient for your question, please ask again.
[lecture #13] 2015.12.10 Control system design for practical system (2/3)†
- return of mini report #2
- review of the experimental apparatus and frequency response experiment
- design example
- determination of plant model(nominal plant and additive uncertainty weight)
nominal.m
subspace.m ... replacement of n4sid in System Identification Toolbox (not provided in IPC)
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
compare_result.m
- room 157 @ Dept. Mech. Bldg.2
ex1203_1
ctrlpref
ex1203_1
346/3.6
ex1203_1
nominal
weight
cont
nominal
compare
compare_result
#ref(): File not found: "2015.12.10-1.jpg" at page "授業/制御工学特論2016"
#ref(): File not found: "2015.12.10-2.jpg" at page "授業/制御工学特論2016"
- Q: Each of the files I wanted to clarify what they represent.
- A: A detailed explanation will be given in the next lecture.
[lecture #14] 2015.12.17 Control system design for practical system (3/3)†
- final report
- design your controller(s) so that the system performance is improved compared with the design example introduced in the previous lecture
- Draw the following figures and explain the difference between two control systems (your controller and the design example):
- bode diagram of controllers
- gain characteristic of closed-loop system from w to z
- time response and frequency spectrum (PSD) of control experiment
- Why is the performance of your system improved(or unfortunately deteriorated)?
- how to improve the performance ?
- order of the nominal plant
- weighting for robust stability
- detailed explanation of m-files in the previous lecture
- specifications of the experimental system
- experimental equipments
- loudspeakers: AURA SOUND NSW2-326-8A (2inch, 15W)
- pressure sensors: NAGANO KEIKI KP15
- A/D, D/A converters: CONTEC AD12-16(PCI), DA12-4(PCI)
- PC: Dell Dimension 1100
- OS: Linux kernel 2.4.22 / Real Time Linux 3.2-pre3
- program sources for frequency response experiment
- freqresp.h
- freqresp_module.c
- freqresp_app.c
- format of spk1.dat (u is used instead of w for spk2.dat)
- 1st column ... frequency (Hz)
- 2nd column ... gain from w(V) to y(V) (signal's unit is voltage (V))
- 3rd column ... phase from w to y
- 4th column ... gain from w to z
- 5th column ... phase from w to z
- program sources for control experiment
- hinf.h
- hinf_module.c
- hinf_app.c
- format of result.dat
- 1st column: time (s)
- 2nd column: z (V)
- 3rd column: y (V)
- 4th column: u (V)
- 5th column: w (V)
- configuration of control experiment
- disturbance signal w is specified as described in hinf.h and hinf_module.c:
#define AMP 3.0 // amplitude for disturbance
#define DIST_INTERVAL 5 // interval step for updating w
count_dist++;
if(count_dist >= DIST_INTERVAL){
w = AMP * (2. * rand() / (RAND_MAX + 1.) - 1.); // uniform random number in [-AMP, AMP]
count_dist = 0;
}
da_conv(V_OFFSET + w, 0); // D/A output to noise source
w is updated with 1ms period (sampling period 0.2ms times DIST_INTERVAL 5)
- control signal u is limited to [-4, 4] as specified in hinf.h and hinf_module.c:
#define U_MAX 4.00
if(u > U_MAX) u = U_MAX;
if(u < -U_MAX) u = -U_MAX;
u is set to 0 for t < 10(s). (controller is operated for 10 <= t < 15.)
- a high pass filter with cut-off frequency are used to cut DC components in z and y as described in hinf.h and hinf_module.c
// HPF(1 rad/s) to cut DC in z and y
#define AF 9.9980001999866674e-01
#define BF 1.9998000133326669e-04
#define CF -1.0000000000000000e+00
#define DF 1.0000000000000000e+00
ad_conv(&yz); // A/D input
// HPFs
yf = CF*xf_y + DF*yz[0];
xf_y = AF*xf_y + BF*yz[0];
zf = CF*xf_z + DF*yz[1];
xf_z = AF*xf_z + BF*yz[1];
experiment directory
- mini exam #2
%-- 12/17/2015 1:36 PM --%
help bodemag
#ref(): File not found: "2015.12.17-1.jpg" at page "授業/制御工学特論2016"
- Q: I wanted to make a print on the final report.
- A: Do you mean you want to submit your report in hard copy? If so, I will receive it.
[lecture #15] 2015.12.24 Control system design for practical system (cont.)†
- preparation of your own controller(s)
- submission procedure of controllers has been changed
- wav file is available
- questionnaires
- to university
- for web-based experimental environment
%-- 12/24/2015 1:11 PM --%
compare
#ref(): File not found: "2015.12.24-1.jpg" at page "授業/制御工学特論2016"