システム制御最適化特論

システム制御最適化特論

担当：平田 健太郎

前期後半 月

限

：

：

号館 第

講義室

7/8

第４回 最適制御

6/17

第１回 最適化問題と線形計画法（

）

第２回 内点法

7/1

第３回 最短経路問題と動的計画法（

）

第４回 最適制御

7/18*

第５回 二次計画法

とモデル予測制御

7/22

第６回 凸解析と線形行列不等式

7/29

第７回 線形行列不等式

による制御系解析・設計

第８回 非線形最適化

* irregular

講義日程（予定）

1

2

3

4

5 50

80

20 10

30

15 15

実行可能基底解が退化

（

解が

より多い）

⇒注意が必要

先の例では

移動は右または上へ

という規則があったため

周回路は存在せず

最短経路が確定したノードの集合を逐次 求めていることが容易であった

これをどう解決するか？

最短経路が確定したノードの集合を



記号の準備

グラフ

は節点集合

と枝集合

からなる

𝑉 = 𝑠, 2,3, ⋯ , 𝑡 , 𝐸 = 𝑖, 𝑗 |

節点

と

の間に枝があるとき



ダイクストラ法

ijkstra Method

多項式時間アルゴリズム

内点法

1

2 4

5 50

15

20 10

30

例題

以下の例に

法を適用してみる

ダイクストラ法の大規模な実行例：

1

2

3

4

5 50

15

80

20 10

30

15

1

2

3

4

5 50

15

80

20 10

30

15

15

1

2

3

4

5 50

80

20 10

30

15

10

30 15

1

2

3

4

5 50

80

20 10

15

30 15

1

2

3

4

5 50

80

20

15

clear;

E=zeros(5,5);

E(1,2)=50;

E(1,3)=80;

E(2,3)=20;

E(2,4)=15;

E(3,4)=10;

E(3,5)=15;

% initialize

S=[];

bS=[1,2,3,4,5];

d=9999*ones(1,5);

d(1)=0;

p=[];

while length(S)<5, S

bS d

nbs=length(bS);

mind=d(bS(1));

v=bS(1);

% Find v

for i=2:nbs,

if d(bS(i))<mind, mind=d(bS(i));

v=bS(i);

end end

5ノードの toy problem に対する Dijkstra法の program code (Octave)

S=[S v];

xbS=[];

for i=1:nbs, if bS(i)~=v,

xbS=[xbS bS(i)];

end end

bS=xbS;

nbs=length(bS);

for k=1:nbs, j=bS(k);

if E(v,j)>0,

if d(j) > d(v)+E(v,j);

d(j)=d(v)+E(v,j);

p(j)=v;

end end end

p

route=5;

prev=p(5);

while prev~=1,

route=[route prev];

prev=p(prev);

end

route=[route prev]

𝑆 の要素を1つずつ増やす： 𝑛 反復高々 n 個の 𝑑 𝑖 を比較

⇒ 𝑂 𝑛^{2 オーダー（最高次の次}_{数のみ係数は無視）}

Dijkstra法の計算量

𝑆 から外に出るエッジについて 𝑑 𝑖 の更新 高々𝑚 回の処理 ⇒ 𝑂 𝑚

(ループ内の処理に見えるが)

エッジの本数は最大でも 𝑚 = 𝑛 𝑛 − 1 2

(𝑛個のノードをすべて相互につなぐとき）

初期化

終了判定 ノード v を探す

(n 個の比較)

n 回のループ 𝑑 𝑗 の付け替え

𝑂 𝑛

𝑂 𝑚

𝑂 𝑚 > 𝑂 𝑛 だから 全体では 𝑂 𝑚𝑛 ?

毎回処理するわけではないので, 𝑂 𝑚𝑛 の見積もりは正しくない. 都合何回あるかを考える.

1

2

3 4

5

1

2

3

4

5

1

2

3

4

5

2 4

∴ 𝑑 𝑗 の付け替えは各ノードに 入ってくる矢印の数だけ行われる.

（行われない場合もある.)

最大でも辺の総数だけ

→

動的計画法（

）

𝑙 𝑖, 𝑗 : 𝑖

から

への最短距離

から直接行ける節点の集合

𝑣∈𝑉_𝑠 𝑙 𝑣, 𝑡 + 𝑎_𝑠𝑣

𝑠

𝑡 𝑣

𝑎_𝑠𝑣

𝑙 𝑣, 𝑡 , 𝑣 ≠ 𝑠 を解くことが 𝑙 𝑠, 𝑡 を解くよりも小さな部分問題に なっていれば, 再帰的に解くことができる(?).

最適性の原理

𝑣∈𝑉_𝑠 𝑎_𝑠𝑣 + 𝑙 𝑣, 𝑡

𝑠

𝑡 𝑣

𝑎_𝑠𝑣

今回の選択

𝑣

までの移動コスト

行った先からゴール までの最小コスト

（最短経路長）

・ 最適制御



最適制御問題



最適制御問題

一般形

非線形）

状態方程式

終端条件の重み関数 状態，制御入力の重み関数 システム

（内部状態 𝑥 𝑡 )

制御入力𝑢(𝑡) 出力

システムの動特性は（非線形な）常微分方程式で与えられるとする

ሶ𝑥 𝑡 = 𝐴𝑥 𝑡 + 𝐵𝑢(𝑡)

線形

𝐽 = 𝑥^𝑇 𝑇 𝑄_𝑓𝑥 𝑇 + න

0 𝑇

𝑥^𝑇 𝑡 𝑄𝑥 𝑡 + 𝑢^𝑇 𝑡 𝑅𝑢 𝑡 𝑑𝑡

一般形

状態方程式

評価関数

連続時間系

線形

(2次評価関数)

一般形

状態方程式

評価関数

𝑥_𝑘+1 = 𝑓 𝑥_𝑘, 𝑢_𝑘 𝑥_𝑘+1 = 𝐴𝑥_𝑘 + 𝐵𝑢_𝑘

離散時間系

𝐽 = 𝐿_𝑓 𝑥_𝑁+1 + ෍

𝑘=0 𝑁

𝐿(𝑥_𝑘, 𝑢_𝑘)

𝐽 = 𝑥_𝑁+1^𝑇 𝑃_𝑁+1𝑥_𝑁+1 + ෍

𝑘=0 𝑁

𝑥_𝑘^𝑇𝑄𝑥_𝑘 + 𝑢_𝑘^𝑇𝑅𝑢_𝑘

動的計画法（

）

離散過程に対する多段階決定問題

状態

に決定

を与えたときの状態遷移：

各段階におけるコスト と終端コストの和

𝑘=1 𝑁

𝐿 𝑥_𝑘, 𝑢_𝑘

を最小にするよう

を選ぶ

𝑥

+𝐿_𝑓 𝑥_𝑁+1

𝑢

初期値 𝑥₁ から始めたとき 𝑁 段決定過程のコスト関数を

𝐽_𝑁^∗ 𝑥₁ で表す. これは初期値の関数なので, 最適コスト関数とよぶ.

𝐽_𝑁 𝑥₁ = ෍

𝑘=1 𝑁

𝐿 𝑥_𝑘, 𝑢_𝑘 + 𝐿_𝑓 𝑥_𝑁+1

とし, 𝐽_𝑁 𝑥₁ を最小化する系列 𝑢_𝑘, 𝑘 = 1, ⋯ , 𝑁 を与えたときのコストを

終端時刻から逆向きに解いていくことを考えると次式が成り立つ

𝐽_𝑁^∗ 𝑥₁ = min

𝑢₁ 𝐿 𝑥₁, 𝑢₁ + 𝐽_𝑁−1^∗ 𝑥₂

𝑥₁ 𝑥₂ 𝑥

𝑘 𝑢

𝑢₁ 𝑘

⋯ ⋯⋯

⋯⋯⋯

𝑥_𝑁+1

𝑢_𝑁

を初期値とする 段決定過程 𝑥₁ 𝑥₂

𝑥

𝑘 𝑢

𝑢₁ 𝑘

⋯ ⋯⋯

⋯⋯⋯

𝑥_𝑁+1

𝑢_𝑁

(最適性の原理)

𝐽_{𝑁−𝑘+1} 𝑥_𝑘 ≔ ෍

𝑖=𝑘 𝑁

𝑥_𝑖^𝑇𝑄𝑥_𝑖 + 𝑢_𝑖^𝑇𝑅𝑢_𝑖 + 𝑥_𝑁+1^𝑇 𝑃_𝑁+1𝑥_𝑁+1

最適コスト関数の初期値 𝐽₀^∗ 𝑥_𝑁+1 = 𝑥_𝑁+1^𝑇 𝑃_𝑁+1𝑥_𝑁+1

𝐽_{𝑁−𝑘+1}^∗ 𝑥_𝑘 = min

𝑢_𝑘,𝑢_𝑘+1,⋯,𝑢_𝑁 𝐽_{𝑁−𝑘+1} 𝑥_𝑘

= 𝑥_𝑘^𝑇𝑄𝑥_𝑘 + 𝑢_𝑘^𝑇𝑅𝑢_𝑘 + 𝐽_𝑁−𝑚 𝑥_𝑚+1

= min

𝑢_𝑘,𝑢_𝑘+1,⋯,𝑢_𝑁 𝑥_𝑘^𝑇𝑄𝑥_𝑘 + 𝑢_𝑘^𝑇𝑅𝑢_𝑘 + 𝐽_𝑁−𝑘 𝑥_𝑘+1

= min

𝑢_𝑘 𝑥_𝑘^𝑇𝑄𝑥_𝑘 + 𝑢_𝑘^𝑇𝑅𝑢_𝑘 +min

𝑢_𝑘+1,⋯,𝑢_𝑁 𝐽_𝑁−𝑘 𝑥_𝑘+1

= min

𝑢_𝑘 𝐿 𝑥_𝑘, 𝑢_𝑘 + 𝐽_𝑁−𝑘^∗ 𝑥_𝑘+1

始めに以下の1段の決定問題を考える.

𝐽₁^∗ 𝑥_𝑁 = min

𝑢_𝑁 𝐿 𝑥_𝑁, 𝑢_𝑁 + 𝐽₀^∗ 𝑥_𝑁+1

= min

𝑢_𝑁 𝑥_𝑁^𝑇𝑄𝑥_𝑁 + 𝑢_𝑁^𝑇𝑅𝑢_𝑁 + 𝑥_𝑁+1^𝑇 𝑃_𝑁+1𝑥_𝑁+1

𝑥_𝑁+1 = 𝐴𝑥_𝑁 + 𝐵𝑢_𝑁

の制約のもとで上記の最小値を与える

は

𝑥_𝑁^𝑇𝑄𝑥_𝑁 + 𝑢_𝑁^𝑇𝑅𝑢_𝑁 + 𝑥_𝑁+1^𝑇 𝑃_𝑁+1𝑥_𝑁+1 = 𝑣^𝑇 𝑅 + 𝐵^𝑇𝑃_𝑁+1𝐵 𝑣 + 𝑥_𝑁^𝑇𝑃_𝑁𝑥_𝑁 𝑣 ≔ 𝑢_𝑁 + 𝑅 + 𝐵^𝑇𝑃_𝑁+1𝐵 ⁻¹𝐵^𝑇𝑃_𝑁+1𝐴𝑥_𝑁

𝑃_𝑁: = 𝐴^𝑇𝑃_𝑁+1𝐴 + 𝑄 − 𝐴^𝑇𝑃_𝑁+1𝐵 𝑅 + 𝐵^𝑇𝑃_𝑁+1𝐵 ⁻¹𝐵^𝑇𝑃_𝑁+1𝐴

𝑢_𝑁 = − 𝑅 + 𝐵^𝑇𝑃_𝑁+1𝐵 ⁻¹𝐵^𝑇𝑃_𝑁+1𝐴𝑥_𝑁 ^{が最小化入力であり} 𝐽₁^∗ 𝑥_𝑁 = 𝑥_𝑁^𝑇𝑃_𝑁𝑥_𝑁

となる

𝐽₂^∗ 𝑥_𝑁−1 = min

𝑢_𝑁−1 𝐿 𝑥_𝑁−1, 𝑢_𝑁−1 + 𝐽₁^∗ 𝑥_𝑁

次のステップは

の最小化であるが, これは先の問題の添字 𝑁 を 𝑁 − 1 にしたもの になっているので, 𝑢_⋅, 𝑃_⋅ を漸化的に更新することにより, 最適解の系列

𝑢₁, 𝑢₂, ⋯ , 𝑢_𝑁 が求められる.

𝑥_𝑘+1 = 𝐴𝑥_𝑘 + 𝐵𝑢_𝑘 𝐽 = ෍

𝑘=0

∞

𝑥_𝑘^𝑇𝑄𝑥_𝑘 + 𝑢_𝑘^𝑇𝑅𝑢_𝑘

無限区間 (Infinite Horizon) の最適制御問題

𝑃 = 𝐴^𝑇𝑃𝐴 + 𝑄 − 𝑃𝐵 𝑅 + 𝐵^𝑇𝑃𝐵 ⁻¹𝐵^𝑇𝑃𝐴 𝑢_𝑘 = − 𝑅 + 𝑃𝐵 ⁻¹𝐵^𝑇𝑃𝐴𝑥_𝑘

Riccati方程式

最適フィードバック則 𝑁 → ∞ のとき, 𝑃_−𝑁 → 𝑃 なる定常解となり

上記の問題は有限区間 (Finite Horizon) を考慮

積分区間を分割

最適性の原理から

動的計画法に基づく連続時間系に対する最適条件の導出

最適コスト関数

𝐽_𝑁 𝑥₁ = min

𝑢₁ 𝐽_𝑁−1 𝑥₂ + 𝐿 𝑥₁, 𝑢₁

最適コスト関数

（離散時間の場合の次式に相当）

𝑥₁ 𝑥₂ 𝑥(𝜏) 𝑥(𝜏₁)

𝜕𝐽^∗

𝜕𝜏 = − min

𝑢 𝜏 ∈𝑈

𝜕𝐽^∗

𝜕𝑥

𝑇

𝑓 𝑥 𝜏 , 𝑢 𝜏 + 𝐿(𝑥 𝜏 , 𝑢 𝜏 )

とし, 右辺第1項をテーラー展開する. （2変数関数）

第2項については

中の を と近似し, の極限をとる.

右辺を最小化する

は

と

𝜕𝑥

に依存するので

これを

𝑢 𝑥 𝜏 ,^𝜕𝐽∗

𝜕𝑥

とおくと

𝜕𝐽^∗

𝜕𝜏 = − min

𝑢 𝜏 ∈𝑈

𝜕𝐽^∗

𝜕𝑥

𝑇

𝑓 𝑥 𝜏 , 𝑢 𝜏 + 𝐿(𝑥 𝜏 , 𝑢 𝜏 )

一般に, 非線形偏微分方程式であるHJB方程式を解析的に解くことは困難

⇒ 制御対象を線形システム, 評価関数を2次形式に限定

さらに最適コスト関数も 𝐽^∗ 𝑥, 𝜏 = 𝑥^𝑇𝑃 𝜏 𝑥 という2次形式に限定する

ሶ𝑥 𝑡 = 𝑓 𝑥, 𝑢 = 𝐴𝑥 + 𝐵𝑢 𝐽 = 𝑥^𝑇 𝑇 𝑄_𝑓𝑥 𝑇 + න

0 𝑇

𝐿 𝑥, 𝑢 𝑑𝑡, 𝐿 𝑥, 𝑢 = 𝑥^𝑇𝑄𝑥 + 𝑢^𝑇𝑅𝑢

下線部のベクトルを 𝑝 とおくと, 括弧内は 𝑝^𝑇 𝐴𝑥 + 𝐵𝑢 + 𝑥^𝑇𝑄𝑥 + 𝑢^𝑇𝑅𝑢

となるので, 平方完成すると

を得る. したがって最小化を達成する 𝑢 （ 𝑢 𝑥, 𝑝ത ） は

𝑝^𝑇 𝐴𝑥 + 𝐵𝑢 + 𝑥^𝑇𝑄𝑥 + 𝑢^𝑇𝑅𝑢 = 𝑢 + 1

2𝑅⁻¹𝐵^𝑇𝑝

𝑇

𝑅 𝑢 + 1

2𝑅⁻¹𝐵^𝑇𝑝 + ⋯

ത

𝑢 𝑥, 𝑝

＝

2𝑅⁻¹𝐵^𝑇𝑝

となる.

𝐽^∗ 𝑥, 𝜏 = 𝑥^𝑇𝑃 𝜏 𝑥 𝜕𝐽^∗

𝜕𝜏 = 𝑥^𝑇 𝑑

𝑑𝜏 𝑃 𝜏 𝑥, 𝜕𝐽^∗

𝜕𝑥 = 𝑝 = 2𝑃𝑥

以上から, 最小化する 𝑢 を代入したあとのHJB方程式は

𝑥^𝑇𝑃𝑥 = −2𝑥ሶ ^𝑇𝑃 𝐴 − 𝐵𝑅⁻¹𝐵^𝑇𝑃 𝑥 − 𝑥^𝑇𝑄𝑥 − −𝑅⁻¹𝐵^𝑇𝑃𝑥 ^𝑇𝑅(−𝑅⁻¹𝐵^𝑇𝑃𝑥)

となる. これが任意の 𝑥 に対して成立するとき

𝑥^𝑇𝑃𝑥 = −2𝑥ሶ ^𝑇𝑃 𝐴 − 𝐵𝑅⁻¹𝐵^𝑇𝑃 𝑥 − 𝑥^𝑇𝑄𝑥 − −𝑅⁻¹𝐵^𝑇𝑃𝑥 ^𝑇𝑅(−𝑅⁻¹𝐵^𝑇𝑃𝑥)

− ሶ𝑃 = 𝑃𝐴 + 𝐴^𝑇𝑃 − 2𝑃𝐵𝑅⁻¹𝐵^𝑇𝑃 + 𝑄 + 𝑃𝐵𝑅⁻¹𝐵^𝑇𝑃

= 𝑃𝐴 + 𝐴^𝑇𝑃 − 𝐵𝑅⁻¹𝐵^𝑇𝑃 + 𝑄

ただし, 両辺はスカラー量なので 2𝑥^𝑇𝑃 𝐴𝑥 = 𝑥^𝑇𝑃𝐴𝑥 + 𝑥^𝑇𝑃𝐴𝑥 ^𝑇

となることを用いた.

よって連続時間線形システムに対する最適制御問題の解は, Riccati微分方程式

− ሶ𝑃(𝑡) = 𝑃(𝑡)𝐴 + 𝐴^𝑇𝑃(𝑡) − 𝑃 𝑡 𝐵𝑅⁻¹𝐵^𝑇𝑃(𝑡) + 𝑄

なる無限区間の最適制御問題となり

その解は

微分方程式 の平衡解

𝑃𝐴 + 𝐴^𝑇𝑃 − 𝑃𝐵𝑅⁻¹𝐵^𝑇𝑃 + 𝑄 = 0, 𝑃 > 0

すなわち代数

方程式の解によって

𝑇 → ∞

とすると

考察の対象は

0

∞

𝑥^𝑇𝑄𝑥 + 𝑢^𝑇𝑅𝑢 𝑑𝑡

なる定数状態フィードバックとして与えられる

＝