授業/制御工学特論2017

授業

Advanced Automation†

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

• outline of this lecture
  • syllabus([https://vos-lc-web01.nagaokaut.ac.jp/])
  • evaluation
    • mini report #1 ... 10%
    • mini exam #1 ... 10%
    • mini report #2 ... 10%
    • mini exam #2 ... 10%
    • final report ... 60%
  • schedule2017 (tentative)
  • map map_v1.1.pdf for review map_v1.1_review.pdf

• review : stabilization of SISO unstable plant by classical and modern control theory
  • transfer functions / differential equations
  • poles / eigenvalues
  • impulse response / initial value response
  • ...

%-- 2017/09/07 13:29 --%
s = tf('s')
P = 1/(s-1)
pole(P)
impulse(P)
Tyr = feedback(P*k, 1)
k = 2
Tyr = feedback(P*k, 1)
P
Tyr = feedback(P*k, 1)
step(Tyr)
k = 0.5
Tyr = feedback(P*k, 1)
step(Tyr)
k = 10
Tyr = feedback(P*k, 1)
step(Tyr)

whiteboard #1

whiteboard #2

whiteboard #3

whiteboard #4

• Q: I don't know how to use the MATLAB.
• Q: マトラボの使い方が分かりません。
• A: A brief explanation of MATLAB and Simulink will be given in the next lecture by using simple examples.
• A: 次週、MATLAB（とSimulink）の使い方について説明します。

• Q: 個人的には倒立振子のようなモデルから制御系の説明が最初にあると制御工学のイメージがつきやすい
• A: 確かに、不安定＝倒れる、という分かり易さは魅力なのですが、1次系で表すことができません。今回は数式の取扱いの容易さから1次系を選びました。 倒立振子については、次回、別の演習で用いていたシミュレーションモデルを紹介したいと思います。

• Q: when you conduct please speak slowly.
• A: Thank you for the sugestion. I will improve my speaking.

[lecture #2] 2017.9.14 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
   • example: stabilization of inverted pendulum (sorry in Japanese)
     • derivation of equation of motion
     • stabilization of 1-link pendulum
     • stabilization of 2-link pendulum

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 ex0914_1.m mod0914_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 ex0914_2.m mod0914_2.mdl
     • Gain Margin(GM); Phase Margin(PM)

%-- 2017/09/14 13:05 --%
a = 1
who
a + 2
demo
lookfor demo
demo toolbox
demo toolbox]
demo toolbox
t = [1, 2, 3]
t = [1 2 3]
u = [1; 2; 3]
t
t'
penddemo
help penddemo
penddemo
ls
ex0914_1
who
P

whiteboard #1

whiteboard #2

whiteboard #3

whiteboard #4

whiteboard #5

• Q: ○7について、 \[ \frac{b_1 s + b_0}{s^2 + a_1 s + a_0} + d \] とした時、対応するSSRは、 \[ \left[\begin{array}{cc|c} 0 & 1 & 0 \\ -a_0 & -a_1 & 1 \\ \hline b_0 & b_1 & 0 \end{array}\right] \] となるという理解で良いでしょうか。（ホワイトボードの説明では、 \[ \left[\begin{array}{cc|c} 0 & 1 & 0 \\ -a_1 & -a_2 & 1 \\ \hline b_1 & b_2 & 0 \end{array}\right] \] に見える為）
• A: その理解で正しいです。赤枠で囲んだ際に添え字が一つずれていることに気付いていませんでした。申し訳ありません。

• Q: I didn't understand the usefulness of SSR.
• A: LQR problem is formulated and solved by using SSR, which will be explain in the next lecture.

• Q: TF から SSR への変換の過程がよくわからなかった。
• A: 与えられた TF に対応する SSR は、状態ベクトルの取り方によって一意に決まらない（無限に存在する）ため、代表的な SSR の型を利用します。 今日扱ったのは controllable canonical form （可制御正準形）で、伝達関数の分母・分子多項式及び直達項をA,C,D行列の特定の場所に当てはめることで、TFからSSRの変換を行ったことになります。この双対の可観測正準形も同様です。他に、今日は紹介しませんでしたが、共振モード毎のブロック対角構造がA行列に表れる型などもあり、目的に応じて使い分けられます。

[lecture #3] 2017.9.21 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
   • ARE and quadratic equation
   • closed loop stability ... Lyapunov criterion
   • Jmin lqr.pdf ≒ proof4.pdf (from B3「動的システムの解析と制御」)
3. example ex0921_1.m mod0921_1.mdl

%-- 2017/09/21 13:17 --%
ex0921_1
B
A
A*B
P
eig(P)
F
J
x0
x0'*P*x0

whiteboard #1

whiteboard #2

whiteboard #3

• Q: MATLABの使い方が所々わからない。
• A: 今日のように時間の余ったときにでも聞いてください。

• Q: QとRの値は何で決まるのですか？
• A: 設計者が決めます。大雑把に言うと、Qを大きく（Rを小さく）すると、状態xの収束性が上がるが大きな制御入力が必要となり、Qを小さく（Rを大きく）すると、状態xの収束性が犠牲になるが小さな制御入力で済む、というトレードオフがあるので、それを考慮して決めます。

• Q: f = b/r P と出来る理由が分からなくなってしまった
• Q: ホワイトボード○5のfとPの関係がどこから出てきたのか分からなかった
• A: LQR問題の解において、状態フィードバックゲインFは代数リカッチ方程式の解Pを用いて与えられます。この関係式（スカラの例題の場合は f = b/r p）を使うと、（リカッチ方程式とは無関係に導出した）f に関する二次方程式が、Pに関する代数リカッチ方程式と同一であることが示されます。わからなければまた聞いてください。

• Q: In m-file, line 13 F = R\B'"P what is this (Yen mark) operator? the same process this? inv(R)\B'

• A: Yen mark means `\'. X\Y and Y/X respectively mean \[ X^{-1} Y, \quad Y X^{-1}. \]

• Q: LQR：L→ 非線形は扱えない？その影響は偏差？時間？パラメータ変動になる？
• A: 先週紹介した状態空間実現は、状態ベクトルに関する一階の線形微分方程式です。よって、非線形特性を直接表現することはできず、平衡点回りで線形化するのが典型的な対処法です。一方、安定となるよう設計された制御系は通常、余裕を持ち、若干の非線形特性は許容されます。この意味で大きな安定余裕を持つように制御系を設計すれば、間接的に非線形特性を考慮できる場合もあります。

• Q: If it is possible, I want to practice more during the lecture, may I get more additional example materials to try on Matlab, for the each topics.
• A: I think this is because a lack of detailed explanation on examples: The mass-spring-damper system given in the today's lecture was unstable one choosing the damping coefficient as a negative value. Please try to change those values; I completely forgot to explain physical meaning of Q and R and how to choose them. As in the Q&A above (sorry in Japanese), matrices Q and R are given by designer considering trade-off between performance and cost e.g. energy consumption. Please try with other values of Q and R, to confirm the trade-off; Therefore, I would like to recommend you to change numerical values in examples for better understandings. If this answer is not sufficient for you, please ask again.

• Q: 少しスピードを落としてほしい
• A: 日本語の説明を少なくしたせいと思います。やはり同じ説明を日本語でもするようにします。

[lecture #4] 2017.9.28 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, $H_2$ control)} \\ \| 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');
   P1 = 1/(s+1);
   bode(P1);
   norm(P1, 'inf')
   P2 = 1/(s^2 + 0.1*s + 1);
   bode(P2);
   norm(P2, 'inf')

4. definition of H infinity norm (SIMO)
5. solve the problem by hand
6. solve the problem by tool(hinfsyn) ex0928.m

%-- 2017/09/28 13:40 --%
s = tf('s');
P1 = 1/(s+1);
bode(P1);
norm(P1, 'inf')
P1
norm(P1, 'inf')
P2 = 1/(s^2 + 0.1*s + 1);
bode(P2)
grid on
norm(P3, 'inf')
norm(P2, 'inf')
format long e
norm(P2, 'inf')
ex0928
ex0929
ex0928
f
K
dcgain(K)
gopt
1/sqrt(29
1/sqrt(2)
gopt

whiteboard #1

whiteboard #2

whiteboard #3

whiteboard #4

whiteboard #5

• Q: H∞制御がまず、どのような設計で使われているなどの具体例があれば嬉しいと思いました
• A: 話が抽象的過ぎてピンとこないという指摘ですね。ただ、この制御対象ならばこの制御理論を使うべき、というルールは特に無く、設計者がどうしたいか、どんな制御系を良しとするかによって制御理論は選んで使われ（無ければ新たに開発され）ます。その意味で、H∞制御が使われる典型的な場面は、制御対象の平均的なモデルとその不確かさ（バラツキ具合、変動）がおおよそ分かっている場合に、それに対する安定性（と性能も）を保証するコントローラ（ロバストコントローラ）の設計です。古い例では小型自動着陸実験機ALFLEXが挙げられます（大野他、ALFLEXのロバスト飛行制御則設計、計測自動制御学会論文集, 34(12):1905--1912（1998））。開発期間が限られており多くの試行錯誤が許されなかったことなどからH∞制御が採用されています。もう一つ、これも古い例ですが、レインボーブリッジや明石海峡橋など、建設中の主塔が強風で振動するのを抑制するためにH∞制御が採用されています（B.F.Spencer Jr. and Michael K.Sain: Controlling Buildings: A New Frontier in Feedback, IEEE Control Systems, 17(6):19{35, 1997）。他にも自動車のサスペンションの制御など、数多くの応用事例があります。調べてみてください。

• Q: 具体的にどのようなシステムを制御する時にH∞、状態FBを使い分けるのか？
• A: 説明が悪かったと思いますが、H∞制御と状態フィードバックは無関係です。つまり、全状態が利用できる状況が状態フィードバック、状態は直接利用できず観測出力を利用する場合が出力フィードバックです。どちらもH∞制御との組み合わせがあります。今日は、手計算では状態フィードバックコントローラをH∞制御で求め、hinfsynでは出力フィードバックH∞コントローラを求めました。出力フィードバックと状態フィードバックをどう使い分けるかという質問に対しては、コストはかかるがセンサが十分設置できるなら制御対象の全状態を利用して性能の良い制御系となるのが状態フィードバックで、多少性能は犠牲になっても制御できれば良いならコストが比較的安い出力フィードバック、ということは言えます。ただ、無限次元系など、センサをいくら設置しても状態フィードバックは不可能な場合もありますが。

• Q: 授業ではふれませんでしたが sup と max の違いを教えていただけないでしょうか？
• A: sup = 最小上界、または上限、max = 最大値、です。1入出力システム P1, P2 \[ P_1(s) = \frac{1}{s+1}, \quad P_2(s) = \frac{s}{s+1} \] を考えると、P1 については ω=0 のゲインがH∞ノルムを与えるので（今日言ったように）sup でも max でも同じですが、P2 については、ω→∞ の極限がH∞ノルム（1）を与えるので、max はとれません。ω=∞ という点がないためです。例えば、十分大きな正実数（10000000とか）をaとすると、 \[ \max_{\omega\in[0,a]} |P_2(j\omega)| = \max_{\omega\in[0,a]} \left|\frac{j\omega}{j\omega+1}\right| = \max_{\omega\in[0,a]} \frac{\omega}{\sqrt{\omega^2 + 1}} = \frac{a}{\sqrt{a^2 + 1}} < 1 \] となって、maxをとる範囲（0からa）を広げればいくらでも真のH∞ノルム 1 に近い値が得られますが、決して1に一致させることはできません。supなら、できます。

• Q: Explanation of H∞ problem was not detailed enough
• A: The difference to LQR is (ii) the target to be minimized: In H infinity control problem, closed-loop H infinity norm is minimized instead of the cost function in LQR. The difference in their concept will be explained in the next lecture. If this answer is not sufficient for you, please ask again.

• Q: What purpose to derive "Tzw" ?
• A: Because our task was to find optimal controller in terms of closed-loop H infinity norm by hand, the f-dependent closed-loop system, Tzw, was first derived. (Then, differentiate with respect to f ...)

• Q: 波形があると分かりやすい
• A: 今日の講義では波形が出ていなくてわかりにくかった、という指摘かと思って回答します：次回、H∞制御とLQRの違いを時間応答で説明します。それぞれで最適な場合として扱われる時間応答です。1回で終わらせられると内容的には区切りが良く分かり易いと思うのですが、分量的に無理のため、現状のような二部構成としています。

[lecture #5] 2017.10.5 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 mod1005.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

%-- 2017/10/05 13:08 --%
sqrt(2-sqrt(2))
mod1005
ex0928
A
B
C
D
h = 0.01
x0 = 0
f = -1+sqrt(2)
zz
h = 0.0001
zz
f = 1
zz
x0 = 1
zz
h
h = 1
x0
x0 = 0
h
zz
sqrt(zz(end)/ww(end))
f
h = 100
sqrt(zz(end)/ww(end))
f = -1+sqrt(2)
sqrt(zz(end)/ww(end))

whiteboard #1

whiteboard #2

whiteboard #3

• Q: H∞のときはノルムを最小としてやっていたが、Jのときは何を最小にしているか分からなかった（コスト関数の意味）
• A: コスト関数（評価関数）の定義は現代制御理論の復習で行ったように \[ J = \int_0^\infty \left\{ x^T Q x + u^T R u \right\} dt \] です。おおざっぱに言えば、Jを最小化することは、x をより速く0に収束させることと、そのためのコスト（制御入力）を小さく済ませることを意味します。しかし通常、両者は相反する（xを速く0に収束させるためには大きなエネルギー（制御入力）が必要となる、または小さなエネルギーで済ませるとxの収束性は悪くなる）ため、そのトレードオフをQとRの調整でとります。という説明が復習の回で不足しておりすみません。

• Q: 最後の証明でγが∞とはどういう意味になるのか？
• A: 今日示した状態フィードバックH∞制御問題の解は、閉ループ系のH∞ノルムがγ未満、つまり \[ \|T_{zw}\|_\infty < \gamma \] を満足するものです。この条件は、γが大きいほど緩和されます。つまり、もしγが1のときに解（コントローラ）が存在すれば、γが2のときも存在します（同じ解が条件を満たします。1は2より小さいので）。逆に言うと、γが1のときに解が存在しなくても、γが2のときには解が存在する可能性があります。2でだめなら100、1000000、とどんどん解が存在する可能性は増え、γ＝∞はつまり、この不等式条件をおかない、最も簡単な問題となります。

• Q: 何を同じにした条件で等価なのか（ふくまれるか）が理解できなかった。
• A: 状態フィードバックH∞制御問題の特別な場合としてLQR問題が含まれることに関しての質問と思います。簡単に言うと、γを∞に持っていき、閉ループ系のH∞ノルムの制約を無くすると、LQR問題になります。

• Q: PID制御と違ってなかなかイメージがわかないと思った。
• A: 次回からH∞制御の intro を3回に分けて説明します。今回と前回のLQRとの関係の話をむしろ後に回した方が、イメージがわいて良かったかもしれません。また感想を教えてもらえると助かります。

• Q: 声が小さいので、もう少し大きな声にしてほしい
• A: 前回まで問題なかったとすれば、今日小さかったのは喉が痛かったからです。次回までには治します。ごめんなさい。

• Q: Why do we need to know the worst case of the disturbance signal w(t) ?
• A: By knowing the worst case w(t) which is very different from the delta function, at least, one can understand two control problems (LQR and H infinity control) are different.

• Q: 何々理解が難しい
• A: 具体的に聞いてもらえると助かります。

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

1. outline: map_v1.1_mixedsens1.pdf
2. H infinity control problem (general case)
3. reference tracking problem
4. weighting function for sensitivity function
5. design example ex1012_1.m ex1012_2.m
6. the small gain theorem
   • proof: Nyquist stability criterion

%-- 2017/10/12 13:38 --%
ex1012_1
P
pole(P)
ex1012_2

whiteboard #1

whiteboard #2

whiteboard #3

whiteboard #4

• Q: should add legend to Matlab plots
• A: I will modify m-files in the remaining lectures.

• Q: ex1012_2.m で行っているのは、Kを最適化しているという理解で良いですか？
• A: 良いです。γの最小化をしているのですが、その説明が抜けておりすみません。

• Q: 重み関数の使用する意味を理解できなかった。
• A: 重み関数 = 周波数依存の重み関数、と説明すべきでした。 もともと |S(j∞)| = 1 のため、周波数依存の重み関数を使用せずにH∞制御問題を解くと（つまり、全ての周波数帯域で |S(jω)| を下げようとすると）K = 0 が γを最小化する最適解（γ=1）となります。説明はしませんでしたが ex1012_2.m の上の方で

  WS = ss(1);

  とすると重み関数を使用しない場合になります。実行してみてください。 一方、低周波数域で大きな値を持つ重み関数を使用してH∞制御問題を解くと（つまり、低周波数域で |S(jω)| を積極的に下げようとすると）、K = 0 は最適解とならなくなり、低周波数域で性能の良いコントローラが得られます。 わからなければまた聞いてください。

• Q: 授業の流れは良いはずなのに、英語と日本語が交互になることで、理解を遅らせている。特に各言語で言い忘れていたりすると良く分からなくなる。
• A: 日本語の分量が少なくなると日本人学生の理解度が悪くなるため、基本的には同じ説明を日本語と英語で二回行うように心がけています。とは言え、厳密に話すことは決めていないため、御指摘のように説明が一対一対応にはなっていないと思います。しかし、指摘の本質は、そもそも日本語と英語の交互の説明自体が分かり難い、という指摘ではないかと思います。これは正直、全く想定していませんでした。様子をみながら対応したいと思いますが、もともと英語の講義ですので、日本語を少なくし、その分、英語で話す際のスピードを落とすしかないかな、と考えています。

• Q: ロバストコントロールに関してよく知らないため、small gain 定理のよさがよくわからなかった。
• A: 次回、ロバスト安定化コントローラの設計法を説明しますので、それで理解できるのではないかと思います。

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

1. review map_v1.1_mixedsens2.pdf and outline
2. normalized uncertainty Delta
3. uncertainty model
4. how to determine P0 and WT
   • example: frequency response of plant with perturbation ex1019_1.m
   • frequency response based procedure for P0 and WT ex1019_2.m
5. robust stabilization problem and equivalent problem
   • design example and simulation ex1019_3.m mod1019.mdl

%-- 2017/10/19 13:30 --%
ex1019_1
ex1019_2
ex1019_3
K
figure(4)
bode(K)
mod1019
c
c = 0.8
c = 0
c = 2

whiteboard #1

whiteboard #2

whiteboard #3

whiteboard #4

• Q: I could only understand ambiguously as a whole. Especially, in the latter half
• A: In righthand side (4) of whiteboard #2, I forgot to explain that P_0 must include all unstable poles of the actual plant \tilde P because the term (1 + W_T \Delta) can represent only stable variations. Therefore, (s-1) is included in the denominator of P_0(s).
  In lefthand side (5) of whiteboard #3, I should have written `for all \tilde P' in the equivalent conditions, e.g., \[ {}^\forall \tilde P \quad \left| \frac{\tilde P(j\omega)}{P_0(j\omega)} - 1 \right| \leq |W_T(j\omega)| {}^\forall \omega \] as confirmed in the gain characteristic plot (some \tilde Ps are specified by varying parameter c).
  I think the last part of the lecture is simple: the robust stabilization problem is formulated with P_0 and W_T after their determination, and the equivalent problem is given by simple block diagram manipulation.
  If you have any questions, please ask again.

• Q: 外乱オブザーバでマイナーループで安定化させるのと、コントローラによって安定化させる違いは？
• A: マイナーループで安定化させる ... 最終的なコントローラがマイナーループとその外側の部分に分かれていて、それぞれを段階的に設計する場合、ということではないかと思います（違っていたら指摘してください）。H∞制御では、一度に丸ごと設計します。別の設計問題に、二自由度制御系があります。閉ループ系の安定性をフィードバックコントローラで確保した後に、フィードフォワードコントローラで所望の目標値追従性能を達成する設計法です。図が無いと分かりにくいと思いますが、一般化プラントの外乱wに目標値信号rを含めると、H∞制御の枠組みで二自由度制御系の二つのコントローラを一度に設計することができます。

• Q: 英語さっぱりわからない。
• A: 特にわからない部分を質問してもらえると助かります。

• Q: \[ \tilde {\bf P} \neq {\bf P} \] ならば \[ \tilde {\bf P} \geq {\bf P} \mbox{ or } \tilde {\bf P} \leq {\bf P} \] とあったが \[ \tilde {\bf P} > {\bf P} \mbox{ or } \tilde {\bf P} < {\bf P} \] になるのでは？
• A: 御指摘の通りです。重み関数の不等式に等号を付けたくてこうしていたのですが、良く考えたら不要でした。すみません。

• Q: モデルの不確かさの決め方、調べ方がよく分からなかった。
• A: 余裕がなくこの辺の話は完全にスキップしました。すみません。（授業では安定性のみで性能については触れませんでしたが）より良い性能を得るために、なるべく現実のプラントを小さな重み関数で、しかし全てを含むようにモデルを構成することが重要です。運動方程式など物理的なモデルがある場合は、その構造に従って相応しいモデルが選べることがあります。そうでない場合、とりあえず加法的摂動モデルや乗法的摂動モデルを使ってみて性能が良くなった方を採用する、ということよく行われます。 劉 康志 著「線形ロバスト制御」計測自動制御学会 編 などが参考になります。
  また、できるだけモデルの次数を低くすることも重要です。得られるコントローラの次数がその分低くなるからです。例えば、 \[ P_0(1+W_T\Delta) = P_0 + (P_0 W_T) \Delta\] より、重み関数を W_T とする乗法的摂動モデルは、重み関数 W を W = P_0 W_T とする加法的摂動モデルとみなせます。しかし、P_0 の次数分、重み関数の次数が高くなるため、このような（P_0と同様の特性が重み関数に含まれる）場合は乗法的摂動モデルの方が有利です。
  不安定極が変動する場合は、LFTの摂動モデルが利用できます（乗法的摂動モデルでは、1 + W_T \Delta の項が安定の場合しか表現できないため、ノミナルプラントに不安定極を含める必要があり、不安定極の変動を直接扱うことはできません。加法的摂動モデルも同様です）。 実際、2入力2出力のシステムMの状態空間実現を \[ M = \left[\begin{array}{c|c:c} 1 & 0.1 & 1 \\ \hline 1 & 0 & 0 \\ 1 & 0 & 0 \end{array} \right] \] として、上側の入出力をΔで閉じると、その伝達関数が \[ \frac{1}{s-1-0.1\Delta} \] となります。Δが-1〜1に変化する実数と思えば、これは不安定極 s = 1 + 0.1Δ を持つ制御対象となります。

[lecture #8] 2017.10.26 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 ex1026_1.m minimize gamma by hand
4. gamma iteration by bisection method ex1026_2.m
5. nominal performance and robust performance ex1026_3.m
6. introduction of robust performance problem

ex1026_1
K
ex1026_2
gam
ex1026_3

whiteboard #1

whiteboard #2

whiteboard #3

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

• schedule2017 mini report #1 and exam #1

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 = [sqrt(2), sqrt(2); -1i/sqrt(2), 1i/sqrt(2)]
     M'
     eig(M'*M)
     svd(M)

   • mini report #1 report1.pdf report1_fixed.pdf ... You will have a mini exam #1 related to this report
4. proof of ||Delta hat||_inf <= 1
5. design example: ex1102_1.m
   • robust performance is achieved but large gap
   • non structured uncertainty is considered ... the design problem is too conservative

M = [sqrt(2), sqrt(2); -1i/sqrt(2), 1i/sqrt(2)]
M'
eig(M'*M)
svd(M)
ex1102_1

whiteboard #1

whiteboard #2

whiteboard #3

whiteboard #4

• Q: "due date" of today's report is wrong.
• A: Thank you for pointing out the mistake, it has been corrected.

• Q: 特異値を使うメリット何ですか？
• A: 1×1の行列（スカラ）の場合は特異値＝絶対値です。特異値は絶対値の自然な拡張となっています。

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

1. return of mini report #1
2. review
   • robust performance but too conservative

     ex1102_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 ex1109_1.m

%-- 2017/11/09 13:07 --%
ex1102_1
M = [1i, 1i; 0, 1]
[U,S,V] = svd(M)
U'*U
V'*V
ex1109_1

whiteboard #1

whiteboard #2

whiteboard #3

whiteboard #4

• Q: ○3の板書は関数a/(s+b)が a<b となっていれば \[ \sup\bar\sigma \left(\left[\begin{array}{cc}\frac{a}{s+b} & 0 \\ \frac{a'}{s+b'} & 0\end{array}\right]\right) \leq 1 \] ということですね
• A: 二度も誤った説明をして大変申し訳ありませんでしたが、a/(s+b) の伝達関数は、直流ゲイン（ω=0における絶対値）が a/b で、これが1より小さかったとしても、 上記の最大特異値は1以下になるとは限りません。実際、a=a'=√2, b=b'=√3 のとき、 \[ \sup\bar\sigma \left(\left[\begin{array}{cc}\frac{a}{s+b} & 0 \\ \frac{a'}{s+b'} & 0\end{array}\right]\right) = \sqrt{\frac{a^2}{b^2} + \frac{{a'}^2}{{b'}^2}} = \sqrt{\frac{4}{3}} \] となり、最大特異値は1より大きくなります。

• Q: What is the circle?
• Q: 最後のシミュレーションの意味が分からなかった。
• A: In the resultant figure of ex1102_1.m, the horizontal and vertical axis indicate the 1st and 2nd elements of vectors (a and b), respectively. For example, if a is given by \[ a = \left[\begin{array}{c} \frac{1}{2}\\\sqrt{3} \end{array}\right], \] then a blue point is displayed at (1/2, √3). In the m-file, vector a is varied so that the point moves around on the unit disk as shown in blue circle. I'm sorry for missing of the explanation.
• A: 上にも書きましたが、横軸と縦軸の意味すら説明していませんでした。ベクトルaまたはbの、第一要素、第二要素をそれぞれ、横軸、縦軸に取ってプロットしました。a は青、bは赤で表示されており、aは単位円周上の座標を使って与えているので、青点は単位円となっています。それぞれのaに対して、b = M a を計算した結果が赤点で、半径が最大となる方向（実線）と、最小となる方向（破線）があることがわかります。このときの入出力の2ノルムの比が、それぞれ最大特異値と最小特異値です。説明不足ですみませんでした。わからなければまた聞いてください。

• Q: \[ \sup_\omega \bar\sigma \left(\left[ \begin{array}{cc} \frac{1}{s+\sqrt{2}} & 0 \\ \frac{1}{s+2} & 0 \end{array}\right]\right) \] sup ω であるが式にωがない。→sをjωに変えるのはどのタイミングでするべきですか？
• A: 上式の s は jω で置き換えるべきでした。僕のミスです。申し訳ありません。

• Q: \[ \|H\|_\infty <1 \mbox{と} \|H\|_\infty < \gamma \] 1以下とγ以下の話の違って？
• A: 板書○1に書いた、R.S.prob for ~Δ （保守的な安定化問題）では、~Δを取り外した閉ループ系のH∞ノルムは1未満であれば良く、これを最小化する必要はありません。 最小化するのは、^Gの中に入っているγで、これを小さくするほど良い制御性能が得られます。つまり、閉ループ系のH∞ノルムが1未満を満たす条件下で、γを最小化する問題を考えています。分からなければまた聞いてください。

• Q: 特異値分解と固有値の関係は？
• A: 今日最後にやった、2×2の実行列に対して考えてみると良いと思います。特異値は、今日解説した通り、入出力ベクトルに関する2ノルムの比の、最大値及び最小値です。一方、固有値の元々の定義は、 \[ M a = λ a \] を満たすような λ です。つまり、入出力の比、という意味はなく、M をスカラ λ と同一視できるような特別な入力ベクトル a が固有ベクトルで、同一視した結果（スカラ λ）が固有値です。

• Q: グラフの赤いバンド幅は、 \[ \|\Delta\|_\infty \leq 1, \quad \|\Delta_P\|_\infty = 1 \] でなく、 \[ \|\tilde \Delta\|_\infty \] のみの大きい幅で考えたため、出てしまう結果という認識でいいでしょうか？ もしそうなら、周波数特性全体で、幅がでないのは何故ですか？
• A: Δを変化させてプロットした、と説明しましたが、ex1102_1.m 中で実際に変化させているのは

  for c=0.8:0.05:1.2
    P = 1/((s-1)*(c*s + 1));
    S_tilde = 1/(1 + P * K_opt);
    bodemag(S_tilde, 'r');
  end

  より、パラメータ c でした。この場合、低周波数域では cs + 1 がほぼ 1 となり、赤の幅はでません。一方、高周波数域ではプラントは P = 1/(s*cs) に近似し、c に依存して赤の幅が出そうな気がしますが、そもそもゲイン|P(jω)| = 1/(cω^2) が0に近づいていくので、S_tilde が 1 に近付き、やはり幅が出なくなります。分からなければまた聞いてください。

[lecture #11] 2017.11.16 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 ex1116_1.m

   ex1102_1
   gam_opt0 = gam_opt
   K_opt0 = K_opt;
   ex1116_1
   gam_opt

5. mini report #2 report2.pdf
   • schedule2017 (modified by blue color text)
6. mini exam #1

%-- 2017/11/16 13:46 --%
ex1102_1
gam_opt
gam_opt0 = gam_opt
Kopt
K_opt
K_opt0 = K_opt
ex1116_1
gam_opt
d
d_opt

whiteboard #1

whiteboard #2

whiteboard #3

• Q: 計算は楽だが、よく理解できないところもあった。
• A: 具体的に聞いてもらえると助かります。

• Q: レポート2の説明のaとbの決め方がわかりません。
• A: Mの中に絶対値が1より大きな要素が含まれる場合、b の中にその要素が含まれるようにa を選べば良いです（bの2ノルムがその大きな要素の絶対値以上となり、aの2ノルムとの比が1を超えることが示せます。ただし、aの2ノルムが1となるように選んだ場合）。 この方法が良くわからなければ、Mの最大特異値を計算し、それが1より大きいことを示しても良いです。

• Q: I did not know why "W" was used, with "W", can γ be made smaller ?
• A: You are right. Roughly speaking, a smaller γ is possible by using a freedom of W as illustrated with today's example.

[lecture #12] 2017.11.30 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 ex1130_1.m
3. introduction of a practical system: active noise control in duct
   • a standard ANC setup
   • 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
   • difficulty in physical modeling
   • experimental setup in lab.
     • symmetrical location of loudspeakers and sensors
       exp_apparatus1.jpg
       exp_apparatus2.jpg
   • frequency response experiment
     ex1130_2.m
     spk1.dat
     spk2.dat
   • room 157 @ Dept. Mech. Bldg. 2

%-- 2017/11/30 13:07 --%
ex1130_1
gamma0
gamma1
format long e
norm(M1_d, 'inf')
ex1130_2
c0 = 346
L = 1.8
c0/(2*L)

whiteboard #1

whiteboard #2

• Q: \[ \frac{c_0}{2L} \] ってどういう意味ですか？ 何故 \[ \frac{c_0}{L} \] ではなく \[ \frac{c_0}{2L} \] なのか？
• A: 音波が一往復して発振した際の共振周波数という意味です。 つまり、2L = 一往復の長さです。 説明が抜けていました。ごめんなさい。

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

1. mini exam #2
   • schedule2017
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
     spk1.dat
     spk2.dat
     
   • 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
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 and frequency spectrum (PSD) of control experiment
   3. Why is the performance of your system improved(or unfortunately deteriorated)?
   • due date: 5th(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 35
   • submit your controller.dat, controller_order.dat, and controller.mat at this page:participant list2017(download is also possible) not later than 27th(Wed) Dec
   • the system will be started until next lecture
   • You can send up to 5 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
   • weighting for robust stability
6. specifications of the experimental system
   1. 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
   2. 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

   3. 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)

   4. 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 0.16Hz (1 rad/s) is used to remove 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];

%-- 2017/12/07 13:01 --%
nominal
n4sid
which n4sid
weight
cont
K
who
Kd
size(Kd.a)
compare
compare_result

whiteboard #1

whiteboard #2

• Q: matlabには元々subspace関数があるため、subspace.mがカレントフォルダにないと、エラーを見ても修正出来ないのを周知した方が良いと思います。
• A: その認識が全くありませんでした。申し訳ありません。次回周知します。

• Q: HPの色が分からない。読めない文字がある。
• A: 色指定を無視して白黒表示にする方法ではどうでしょうか。または、見にくい色を教えてもらえれば適宜変更します。

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

• return of mini exam #2
• review
  • Q/A: original Matlab function subspace() is hidden by subspace.m (please locate subspace.m in the same working directory for controller design.)
  • 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)
• preparation of your own controller(s) by using the remote experiment system

%-- 2017/12/14 13:47 --%
compare
compare_result_fixed
help print
compare_result_fixed

whiteboard #1

• Q: 1/5に提出するレポートについては、ファイル形式に指定はありすでしょうか。（pdfなど）
• A: pdfまたはdoc形式で提出してください。このページに

  submit your report(pdf or doc) by e-mail to kobayasi@nagaokaut.ac.jp

  と記載している通りです。

• Q: 遠隔実験用の web ページに登録した際の"Registration Finished" のページで Password が全部表示されるのはあまりよくないと思う。（他と同じPWを使う人もいると思うので）
• A: 御指摘ありがとうございます。席が近いのを考慮していませんでした。近日中に対応させていただきます。

• Q: 日本語がないと苦しい
• A: 内容が具体的なのでわかるだろう、とかなり勢いで喋っていたところがあったかもしれません。御指摘に感謝します。

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

• preparation of your own controller(s)
• questionnaires
  • to university
  • for web-based experimental environment

%-- 2017/12/21 13:32 --%
nominal
weight
compare
weight
cont
compare
weight
cont
compare
nominal
weight
cont
compare
weight
compare_result
cont

• Q: 全体として日本語で講義が行われたらより良かった。
• A: ご意見ありがとうございます。そもそも分かり易い英語でなかったと思います。まずはこの点を改善します。

Contents

0．トップ

1. 掲示板

2．メンバー

3. 研究テーマ

4. 研究生活

5．イベント

6．授業

7．学生実験

8．研究室内のページ

9．情報公開

最新の10件