ディジタル制御（1987年）その5 v2 pythonで記述してみる（7.固有値を指定した制御系）

2024.4.7 ディジタル制御（1987年）その5 v2 pythonで記述してみる（7.固有値を指定した制御系） 下 条 誠 電気通信大学名誉教授 https://researchmap.jp/read0072509/ https://www.docswell.com/user/m_shimojo “高橋安人，ディジタル制御 ，岩波書店， 1987 The University of Electro-Communications https://github.com/m4881shimojo/YTakahashi_Digital-control/tree/main Department of Mechanical Engineering and Intelligent System

• https://researchmap.jp/read0072509/
• https://www.docswell.com/user/m_shimojo
• https://github.com/m4881shimojo/YTakahashi_Digital-control/tree/main

目 次 1. 緒 言 1.1 ディジタル制御工学 4.4 4.5 1.2 1.3 1.4 1.5 5. 連続時間プラントの離散時間形 5.1 フィードバック制御系 5.2 １次プラントの制御 ディジタル制御系の構成 量子化、サンプリングおよびホールド 例－ロボットアームのデッドビート応答 プログラム作成入門 2. ｚ変換と反復計算 2.1 ｚ領域における時系列 2.2 入力と出力 2.3 極とゼロ 2.4 2.5 調和振動 リカーシブ・アルゴリズム 5.3 5.4 5.5 時系列からの伝達関数 観測器による状態推定 追従系の設計 PID制御則 PIDゲイン調整 6. 線形2乗最適制御 6.1 行列リカチ(Riccati)式 6.2 慣性体のLQ制御 6.3 I動作を含むLQ制御 6.4 観測器を含むLQI制御 6.5 ０型プラントのLQI制御 ３. 線形離散時間系 3.1 3.2 3.3 状態空間系 伝達関数 同伴形 3.4 応答の算定 7.1 LQ制御系の固有値 3.5 安定に関する定理 7.2 7.3 7.4 7.5 LQ制御系の根軌跡 固有値を指定の制御系および状態推定系設計 有限整定系 多変数制御系 ４. 連続時間プラントの離散時間形 4.1 プラント微分方程式と等価の差分式 4.2 パルス伝達関数と拡張ｚ変換 4.3 解析的手法によるｚ変換と拡張ｚ変換 7. 8. 固有値を指定した制御系 同定および適応制御 2 8.1 8.2 簡単な適応サーボ系 模型規範適応制御 8.3 8.4 8.5 MRASによる同定と適応制御への応用 模型追従適応制御 最小2乗同定アルゴリズム 9. 周波数領域における諸関係 9.1 振動状定常状態 9.2 周波数応答 9.3 制御系の周波数応答 9.4 観測器を含む状態ベクトルフィードバック系 9.5 離散形式のフーリエ変換 10. ランダム信号系 10.1 乱数から造成する正規性白色ノイズ 10.2 10.3 10.4 有色ノイズと自己相関 共分散行列 連続時間ノイズを含む離散時間形 10.5 最適フィルタ

コメント 3 本書の面白いところは、全てのアルゴリズムが、HP-Basi言語と付属の基本行列 演算で作られている点です。正準形変換の方法やリカチ行列の解法などの記述が 多々あり大変興味深いものです。これらは現在ではpython-controlなどを使えば 簡単に解けます。しかし、その解法の基本原理まで踏み込み、意外と簡単なアル ゴリズムで計算できることを本書を通じて実感できることは楽しいです。なお本 解説は本人の再学習のためであり、主な内容のメモとプログラムの記述だけです。 デジタル制御の解説ではありません。 また、作成したpythonプログラムにはミスなどあると思います。書籍の結果とは 比べチェックはしましたが完全ではありません。その点ご注意ください。 （Pythonのprogramは下記のGitHubに格納してあります） https://github.com/m4881shimojo/YTakahashi_Digital-control/tree/main 高橋 安人氏（1912-1996)は、日本の機械工学者。工学博士（東京帝国大学）、 カリフォルニア大学バークレー校名誉教授、豊橋技術科学大学名誉教授。当初は 伝熱工学を専攻していたが、太平洋戦争後は日本の制御工学の草分けとして活躍 した。https://ja.wikipedia.org/wiki/高橋 安人

• https://github.com/m4881shimojo/YTakahashi_Digital-control/tree/main
• https://ja.wikipedia.org/wiki/

7. 7.1 固有値を指定した制御系 LQ制御系の固有値 最適(LQ,LQI)制御系の根軌跡を求める 7.2 LQ制御系の根軌跡 7.3 固有値を指定の制御系および状態推定系設計 7.4 有限整定系 7.5 多変数制御系 指定した通りの閉ループ固有値（極） をもつ制御系を設計することが主題。 最適制御系（LQ,LQI)の根分布が固有値 指定の参考になる 重み係数wを変化させたときの根軌跡を 調べる 4

7.1 LQ制御系の固有値 5 wを変化→最適制御系の固有値の根軌跡は？ ∞ 評価関数： 𝐽 = ෍ 𝑊𝑦 𝑦 2 𝑘 + 1 + 𝑤𝑢2 𝑘 (6-2) 𝑘=0 最適(LQ,LQI)制御系の 重みｗを変数とした根軌 跡を求める wをパラメータとした根軌跡

7.1 LQ制御系の固有値 6 最適(LQ,LQI)制御系の重みｗを変数とした根軌跡 最適制御系（LQ, LQI)の根軌跡を求めるため、補助変数 h を導入する LQ制御系の特性式 補助ベクトルhの導入： 𝒙 𝑘 + 1 = 𝑷𝒙 𝑘 − 𝒒𝑤 −1 𝒒𝑇 𝒉(𝑘) (7-4) 𝒉 𝑘 − 1 = 𝑾𝑥 𝒙 𝑘 + 𝑷𝑇 𝒉(𝑘) 𝒙 ∈ ℝ𝑛×1 , 𝒉 ∈ ℝ𝑛×1 , 𝒒 ∈ ℝ1×𝑛 𝑷 ∈ ℝ𝑛×𝑛 , 𝑾𝑥 ∈ ℝ𝑛×𝑛 上記（連立ベクトル差分式）の特性式 Fc(z)： 𝑧𝑰 − 𝑷 𝐹𝑐 𝑧 = 𝑑𝑒𝑡 −𝑾𝑥 𝒒𝑤 −1 𝒒𝑇 𝑧 −1 𝑰 − 𝑷𝑇 𝑤∈ℝ （システム制御のためのマトリクス理論 = 𝑑𝑒𝑡 𝑧𝑰 − 𝑷 ∙ 𝑑𝑒𝑡 𝑧 −1 𝑰 − 𝑷𝑇 + 𝑾𝑥 𝑧𝑰 − 𝑷 = 𝑑𝑒𝑡 𝑧𝑰 − 𝑷 ∙ 𝑑𝑒𝑡 𝑧 −1 𝑰 − 𝑷𝑇 ∙ 𝑑𝑒𝑡 𝑰 + 𝑾𝑥 𝑧𝑰 − 𝑷 A z 開ループ系特性式 A 𝑧 −1 𝐺 𝑧 = 𝑐 𝑧𝑰 − 𝑷 −1 −1 p41参照） 𝒒𝑤 −1 𝒒𝑇 𝒒𝑤 −1 𝒒𝑇 𝑧 −1 𝑰 − 𝑷𝑇 𝑑𝑒𝑡 1 + 𝑤 −1 𝐺 𝑇 𝑧 −1 𝐺(𝑧) −1 𝒒 ：閉ループ伝達関数 −1 (e)

7.1 LQ制御系の固有値 1入出力系とすると LQ制御系の特性式 (e)式は上記の ようになる 𝐵(𝑧) 𝐺 𝑧 = , 𝐴(𝑧) 𝐹𝑐 𝑧 = 𝐴 𝑧 ∙ 𝐴 𝐺 𝑧 𝑧 −1 −1 7 𝐵(𝑧 −1 ) = 𝐴(𝑧 −1 ) 1 + 𝐵 𝑧 ∙ 𝐵 𝑧 −1 𝑤 (7-5) （１）重みw→∞ならば、u(k)=0（無制御）が最適制御となる 𝐹𝑐 𝑧 = 𝐴 𝑧 ∙ 𝐴 𝑧 −1 (7-6) （２）重みw→0ならば、Eq.7-7漸近しゼロ(伝達関数分子多項式の根)が閉ループ極となる 𝐹𝑐 𝑧 = 𝐵 𝑧 ∙ 𝐵 𝑧 −1 次節では、重みwの値を変えたときの根軌跡を検討する (7-7)

7.2 LQ制御系の根軌跡 7.1 LQ制御系の固有値 最適(LQ,LQI)制御系の根軌跡を求める 7.2 LQ制御系の根軌跡 7.3 固有値を指定の制御系および状態推定系設計 7.4 有限整定系 7.5 多変数制御系 指定した通りの閉ループ固有値（極） をもつ制御系を設計することが主題。 最適制御系（LQ,LQI)の根分布が固有値 指定の参考になる 重み係数wを変化させたときの根軌跡を 調べる 8

7.2 LQ制御系の根軌跡（例1） LQ 制御系での重みwを、パラメータとする根軌跡を示す 例1の対象システム： 𝑦(𝑡) 運動方程式： u(t) 𝑦ሷ 𝑡 + 𝑦(𝑡) ሶ =𝑢 𝑡 b m m=1 b=1 ｚ変換式： （5.3 𝐺𝑝 𝑧 = 追従系の設計、6.2慣性体のLQ制御を参照） 状態方程式： 𝒙 𝑘 + 1 = 𝑷𝒙 𝑘 + 𝒒𝒖 𝑘 𝑷= 1 0 1−𝑝 𝑝 𝑝 = 𝑒 −𝑇 , 𝒒= 𝑝+𝑇 1−𝑝 𝑎 = 𝑝 + 𝑇 − 1, 𝑒𝑥) 𝑇 = 0.1 ∶ 𝑝 = 0.904837, 𝑎𝑧 + 𝑏 𝐵(𝑧) = 𝑧−1 𝑧−𝑝 𝐴(𝑧) 𝒄= 1 𝒚 𝑘 = 𝒄𝒙(𝑘) 0 𝑇: サンプリング周期 𝑏 = 1 − 𝑝 − 𝑝𝑇 𝑎 = 0.004837 𝑏 = 0.004679 9

7.2 LQ制御系の根軌跡（例1） LQ制御系の特性式： 𝐺𝑝 𝑧 = 𝐹𝑐 𝑧 = 𝐴 𝑧 ∙ 𝐴 𝑧 −1 + 𝐹𝑐 𝑧 = 0 10 𝑎𝑧 + 𝑏 𝐵(𝑧) = 𝑧−1 𝑧−𝑝 𝐴(𝑧) 1 𝐵 𝑧 ∙ 𝐵 𝑧 −1 𝑤 𝑧 4 + 𝑎(1)𝑧 3 + 𝑎(2)𝑧 2 + 𝑎 3 𝑧 + 𝑎(4) = 0 𝑎𝑏 − 𝑤 1 + 𝑝 𝑎 1 = 𝑤𝑝 2 𝑤 1 + 𝑝2 + 1 + 𝑝 𝑎 2 = 𝑤𝑝 𝑎𝑏 − 𝑤 1 + 𝑝 𝑎 3 = 𝑤𝑝 𝑎 4 =1 2 2 + 𝑎2 + 𝑏2 (7-9) 𝑎, 𝑏, 𝑝 ← 𝑔𝑖𝑣𝑒𝑛 𝑤 ← 𝑝𝑎𝑟𝑎𝑚𝑒𝑡𝑒𝑟 このwをパラメータとする4次方程式を解き、 wに関する根軌跡を次頁に示す

7.2 LQ制御系の根軌跡（例1） 11 LQ制御系特性式の極とゼロ 𝑎𝑧 + 𝑏 𝐵(𝑧) 𝐺 𝑧 = = 𝑧−1 𝑧−𝑝 𝐴(𝑧) 𝑒𝑥) 𝑇 = 0.1 𝑝 = 0.904837 𝑎 = 0.004837. LQ制御系の特性式： 𝐹𝑐 𝑧 = 𝐴 𝑧 ∙ 𝐴 𝑧 −1 + 分母を払い記述：𝑤 𝑧 2 − 1 + 𝑝 𝑧 + 𝑝 １． w→∞ 𝑏 = 0.004679 （5.3節参照） 1 𝐵 𝑧 ∙ 𝐵 𝑧 −1 𝑤 𝑝𝑧 2 − 1 + 𝑝 𝑧 + 1 + 𝑧(𝑎𝑧 + 𝑏)(𝑏𝑧 + 𝑎) =0 A(z)によって決まる 𝐹𝑐 𝑧 = 𝐴 𝑧 ∙ 𝐴 𝑧 −1 = 0 ２． w→０ (7.-8) 𝐴 𝑧 = 𝑧−1 𝑧−𝑝 𝒛=𝟏, 𝟎.𝟗𝟎𝟒𝟖𝟑𝟕 𝐴 𝑧 −1 = 𝑧 −1 − 1 𝑧 −1 − 𝑝 𝒛=𝟏, 𝟏.𝟏𝟎𝟓𝟏𝟕 B(z)によって決まる 𝐹𝑐 𝑧 = 𝐵 𝑧 ∙ 𝐵 𝑧 −1 = 0 𝑧 𝑎𝑧 + 𝑏 𝑏𝑧 + 𝑎 = 0 -a/b=-1.0337678991237444, -b/a= -0.9673351250775274 𝑧=0, -1.03377, -0.96734

7.2 LQ制御系の根軌跡（例1） : w→∞ : w→0 w=1e-05 12 例1）最適制御系（LQ) の重みｗをパラメータ とする根軌跡 6.2節で設計した慣性体 のLQ制御系 単位円 • 単位円が安定限界 • • 閉ループ極は一対が円内 （安定根） 他の一対が円外（不安定 • 根）にある それらの軌跡はw→0の とき閉ループゼロ（図中 〇）へ漸近する wの値（始まりを示す） DGC_no5/ p125v2

7.2 LQ制御系の根軌跡（例1） : w→∞ : w→0 単位円 13 例1）最適制御系（LQ) の重みｗをパラメータ とする根軌跡 （拡大図: 単位円内） • • 単位円が安定限界 閉ループ極は一対が円内 （安定根） • • 他の一対が円外（不安定 根）にある それらの軌跡はw→0の とき閉ループゼロ（図中 〇）へ漸近する LQ制御では単位円内 の固有値を使う DGC_no5/ p125v2

7.2 LQ制御系の根軌跡（例1） 14 例1）最適制御系（LQ) の重みｗをパラメータ とする根軌跡 : w→∞ : w→0 単位円 （さらに拡大図） w=0.001 w→∞ • w=0.1 • • • w の 値 始 ま り を 示 す (1, 0.904, 1, 1.105) の4点から始まり 単位円内で分岐する対と 単位円外で分岐する対が ある それらの軌跡はw→0の とき閉ループゼロへ漸近 する (1.0にある×点は重根) DGC_no5/ p125v2Asobi

7.2 LQ制御系の根軌跡（例1） page125 “高橋安人,ディジタル制御 ,岩波書店,1987 15

7.2 LQ制御系の根軌跡（例2） 例２の 対象システム： 𝐺𝑝 𝑧 = 6.3節で既出 𝑏1 𝑧 + 𝑏2 𝑧 − 𝑝1 𝑧 − 𝑝2 𝑏1 = 1 − 𝑝1 − 𝑟𝑝2 Τ 1 − 𝑟 𝑝1 = 𝑒 −𝑇Τ𝑇1 − 𝐾0 + − 𝑟 = 𝑇2 Τ𝑇1 ≠ 1 𝐺𝑝 𝑧 V(𝑧) 𝑰: 𝑧 𝑧−1 𝑝2 = 𝑒 −𝑇Τ𝑇2 （4.3節参照） (4-12) 𝑏2 = 𝑝1 𝑝2 − 𝑝2 − 𝑟𝑝1 Τ 1 − 𝑟 𝑅(𝑧) + 𝐸(𝑧) 16 + + 𝑏1 𝑧 + 𝑏2 𝑧 − 𝑝1 𝑧 − 𝑝2 U(𝑧) Y(𝑧) 𝑲 (6.3 LQI制御系) I動作を直列結合すると開ループ伝達関数は次のようになる 𝑧 𝑏1 𝑧 + 𝑏2 𝐵(𝑧) G 𝑧 = = 𝑧 − 1 𝑧 − 𝑝1 𝑧 − 𝑝2 𝐴(𝑧)

7.2 LQ制御系の根軌跡（例2） 例２の 対象システム + I 動作 G 𝑧 = 17 𝑧 𝑏1 𝑧 + 𝑏2 𝐵(𝑧) = 𝑧 − 1 𝑧 − 𝑝1 𝑧 − 𝑝2 𝐴(𝑧) LQI制御系の特性式： (7-5)式を適用、分母分子多項式を整理する 𝐹𝑐 𝑧 = 𝑧 6 + 𝑎(1)𝑧 5 + 𝑎(2)𝑧 4 + 𝑎(3)𝑧 3 + 𝑎(4)𝑧 2 + 𝑎 5 𝑧 + 𝑎(6) = 0 𝑎 1 = −2 − 𝑝1 + 1Τ𝑝1 − 𝑝2 + 1Τ𝑝2 𝑎 2 = 3 + 2 𝑝1 + 1Τ𝑝1 + 2 𝑝2 + 1Τ𝑝2 + 𝑝1 + 1Τ𝑝1 𝑝2 + 1Τ𝑝2 − 𝑏1 𝑏2 Τ 𝑤𝑝1 𝑝2 𝑎 3 = −4 − 2 𝑝1 + 1Τ𝑝1 − 2 𝑝2 + 1Τ𝑝2 − 2 𝑝1 + 1Τ𝑝1 𝑝2 + 1Τ𝑝2 − 𝑏12 + 𝑏22 Τ 𝑤𝑝1 𝑝2 𝑎 4 =𝑎 2 𝑎 5 =𝑎 1 𝑎 6 =1 6次の多項式で安定根はその半分の３つある。

7.2 LQ制御系の根軌跡（例2） LQI制御系特性式の極とゼロ 𝑧 𝑏1 𝑧 + 𝑏2 𝐵(𝑧) G 𝑧 = = 𝑧 − 1 𝑧 − 𝑝1 𝑧 − 𝑝2 𝐴(𝑧) LQI制御系の特性式： １．w→∞ 𝐹𝑐 𝑧 = 𝐴 𝑧 ∙ 𝐴 𝑧 −1 + 1 𝐵 𝑧 ∙ 𝐵 𝑧 −1 𝑤 A(z)によって決まる極 𝐴 𝑧 = 𝑧−1 𝑧 − 𝑝1 𝑧 − 𝑝2 = 0 → 𝑧 = 1, 𝑝1 , 𝑝2 𝐴 𝑧 −1 = 𝑧 −1 − 1 ２．w→０ 𝑒𝑥) 𝑇 = 4 𝑝1 = 0.670320, 𝑝2 = 0.5312417 𝑏1 = 0.094326, 𝑏2 = 0.066091 𝑧 −1 − 𝑝1 𝑧 −1 − 𝑝2 = 0 → 𝑧 = 1, 1/𝑝1 , 1/𝑝2 B(z)によって決まるゼロ B 𝑧 = 𝑧 𝑏1 𝑧 + 𝑏2 = 0 → 𝑧 = 0, - 𝑏2 / 𝑏1 B 𝑧 −1 = 𝑧 𝑏1 𝑧 −1 + 𝑏2 = 0 → 𝑧 = ∞, -𝑏1 / 𝑏2 z = 1.0 𝑧 = 1.0 𝑝1 = 0.670320 𝑝2 = 0.5312417 1Τ𝑝1 = 1.49182 1Τ𝑝2 = 1.88238 𝑧=0 𝑧=∞ −𝑏2 Τ𝑏1 = −0.70067 −𝑏1 Τ𝑏2 = −1.42721 18

7.2 LQ制御系の根軌跡（例2） : w→∞ 19 例２）最適制御系 （LQI)の重みｗをパラ メータとする根軌跡 : w→0 • 6次の多項式では根が６つある • 実軸１では重根である • 安定根は、その半分の３つある 単位円 • 根軌跡はw→∞において開ループ 極×から出発して、 w→0にとも 1.49182 1.8823 −1.42721 ない開ループのゼロ〇へ漸近する。 極：X ゼロ：〇 z = 1.0 𝑧 = 1.0 𝑝1 = 0.670320 𝑝2 = 0.5312417 1Τ𝑝1 = 1.49182 1Τ𝑝2 = 1.88238 𝑧=0 𝑧=∞ −𝑏2 Τ𝑏1 = −0.70067 −𝑏1 Τ𝑏2 = −1.42721 DGC_no5/ p126v1

7.2 LQ制御系の根軌跡（例2） 20 単位円内について : w→∞ : w→0 • 根軌跡はw→∞において開ループ極× から出発する • ｗを小さくするにつれて、B,C極から 出発した軌跡は複素極へと分岐し卵 単位円 型の経路を経て負の実軸で合流し、 w→0にともない開ループのゼロ〇へ 漸近する w=0.01 w=1 A B • 原点では〇が重複する。その一つは C I動作z/(z-1)によるゼロ〇である。 • 一般にn次でゼロの数がmのシステム -1.427 では、(n-m)個の根がｚ＝０へ漸近 -0.70067 する。この例ではn=3,m=2である • このLQI制御では、 w→0のとき閉 ループ極の2つがz=0もう一つがz=0.70067となり、後者が開ループゼ ロと極ゼロ相殺をおこす • 従って閉ループ系の応答はz=0にあ 𝐴 = 0.670320 𝐵 = 0.5312417 𝐶 = 1.0 る2根だけに支配され支配され、有限 整定応答を与える DGC_no5/ p126v1

7.2 LQ制御系の根軌跡（例2） 21 拡大図 • このLQI制御では、 w→0のとき閉 ループ極の2つがz=0、もう一つが 単位円 z=-0.70067となり、後者が開ルー プゼロと極ゼロ相殺をおこす • 従って閉ループ系の応答はz=0に ある2根だけに支配され、有限整定 応答を与える。 • この有限整定は(n-m)刻みを要する。 ただし、LQまたはLQI制御でw→0 にすると必ずしも有限整定になる 0.5312417 0.670320 1.49182 とは限らない。開ループゼロが単 位円外にあるシステムではそうな らない（不安定ゼロ） • なぜならw→0で安定根軌跡の1つ が開ループゼロの逆数へ漸近する ので、開ループゼロとの相殺がお こらないためである • このように開ループゼロは制御系 の挙動に大きな影響を及ぼす。 (書籍p.127) DGC_no5/ p126v1

7.2 LQ制御系の根軌跡（例2） 22 （4.3節参照） page126 “高橋安人,ディジタル制御 ,岩波書店,1987

7.3 固有値を指定の制御系および状態推定系設計 7.1 LQ制御系の固有値 7.2 LQ制御系の根軌跡 好ましい挙動を示す 固有値の見当をつける 7.3 固有値を指定の制御系および状態推定系設計 7.4 有限整定系 7.5 多変数制御系 固有値を指定した制御系 の設計法と状態推定系を 双対問題として取り扱う 23

7.3 固有値を指定の制御系および状態推定系設計 プラントがn次の伝達関数とする 好ましい挙動を示す固有値の見当がついた後の話し 𝑏1 𝑧 𝑛−1 + 𝑏2 𝑧 𝑛−2 + ⋯ + 𝑏𝑛−1 𝑧 + 𝑏𝑛 𝐵(𝑧) 𝐺 𝑧 = 𝑛 = 𝑧 + 𝑎1 𝑧 𝑛−1 + 𝑎2 𝑧 𝑛−2 + ⋯ + 𝑎𝑛−1 𝑧 + 𝑎𝑛 𝐴(𝑧) (7-14) プラントが0型→I(積分)動作つける→全体次数(n+1)となる (n+1)個の固有値を指定すると次の閉ループ特性式となる 閉ループ特性式： 𝐹𝑐 𝑧 = 𝑧 𝑛+1 + 𝛼1 𝑧 𝑛 + 𝛼2 𝑧 𝑛−1 + ⋯ + 𝛼𝑛 𝑧 + 𝛼𝑛+1 (7-15) (7-14)を可制御同伴形（可制御正準形）で表すと 𝒙𝑐 𝑘 + 1 = 𝑷𝑐 𝒙𝑐 𝑘 + 𝒒𝑐 𝑢(𝑘) 0 0 0 𝑷𝒄 = ⋮ 0 −𝑎𝑛 1 0 0 ⋮ 0 −𝑎𝑛−1 0 1 0 ⋮ 0 ⋯ … ⋯ ⋯ ⋮ ⋯ ⋯ 0 0 0 ⋮ 1 −𝑎1 24 0 0 0 𝒒𝒄 = ⋮ 0 1 𝒚 𝑘 = 𝑪𝒙𝑐 𝑘 𝑪 = 𝑏𝑛 𝑏𝑛−1 𝑏𝑛−2 ⋯ 𝑏1

7.3 固有値を指定の制御系および状態推定系設計 25 このプラントへゲイン行列 kc の状態フィードバックをつけると 0 0 0 𝑷𝒄 − 𝒒𝒄 𝒌𝒄 = ⋮ 0 −𝑎𝑛 − 𝑘1 1 0 0 ⋮ 0 −𝑎𝑛−1 − 𝑘2 0 1 0 ⋮ 0 ⋯ … ⋯ ⋯ ⋮ ⋯ ⋯ 0 0 0 ⋮ 1 −𝑎1 − 𝑘𝑛 𝒌𝑐 = 𝑘1 𝑘2 𝑘3 ⋯ 𝑘𝑛 伝達関数は、B(z)は(7-14)と同じだが、分母A(z)は下記となる B(z) 𝐺𝑓 𝑧 = 𝐴𝑓 z = 𝑧 𝑛 + 𝑎1 + 𝑘𝑛 𝑧 𝑛−1 + ⋯ + 𝑎𝑛−1 − 𝑘2 𝑧 + 𝑎𝑛 + 𝑘1 𝐴𝑓 (z) プラントが0型 → I(積分)動作つける必要がある。すると開ループ伝達関数は 𝑘0 𝑧 𝐵(𝑧) G0 𝑧 = 𝑧 − 1 𝐴𝑓 (𝑧) 従って、閉ループ特性式は1+G0(z)の 分母を払って次のように求められる

7.3 固有値を指定の制御系および状態推定系設計 26 従って、閉ループ特性式は1+G0(z)の分母を払って次のように求められる 𝐹𝑐 𝑧 = 𝑧 − 1 𝐴𝑓 𝑧 + 𝑘0 𝑧𝐵 𝑧 = 𝑧 𝑛+1 + 𝑎1 + 𝑘𝑛 − 1 + 𝑘0 𝑏1 𝑧 𝑛 + + 𝑎2 + 𝑘𝑛−1 − 𝑎𝑛 + 𝑘1 + 𝑘0 𝑏2 𝛼2 𝑧 𝑛−1 + ⋯ + 𝑎𝑛 + 𝑘1 − 𝑎𝑛−1 + 𝑘2 + 𝑘0 𝑏𝑛 𝑧 + − 𝑎𝑛 + 𝑘1 (7-16) この特性式と(7-15)式の係数を対応項ごとに比較すると 𝛼1 = 𝑎1 + 𝑘𝑛 −1 + 𝑘0 𝑏1 𝛼2 = 𝑎2 + 𝑘𝑛−1 − 𝑎𝑛 + 𝑘1 + 𝑘0 𝑏2 ⋯⋯⋯⋯ 𝛼𝑛 = 𝑎𝑛 + 𝑘1 − 𝑎𝑛−1 + 𝑘2 + 𝑘0 𝑏𝑛 𝛼𝑛+1 = − 𝑎𝑛 + 𝑘1 (7-17) (7-17)式の両辺の総和を行う と相殺がおこり、k0が判明する。 これを(7-17)式の下から上に 順次代入して求める 𝑘0 = 1 + 𝛼1 + ⋯ + 𝛼𝑛+1 1 + 𝑏1 + ⋯ + 𝑏𝑛 (7-18) 𝑘1 = −𝛼𝑛+1 − 𝑎𝑛 𝑘2 = 𝑎𝑛 + 𝑘1 − 𝛼𝑛 − 𝑎𝑛−1 + 𝑘0 𝑏𝑛 ⋯⋯⋯⋯ 𝑘𝑛 = 𝑎2 + 𝑘𝑛−1 − 𝛼2 − 𝑎1 + 𝑘0 𝑏2 (7-19)

7.3 固有値を指定の制御系設計Ackermann公式 27 プラントが伝達関数(7-14)や可制御同伴形（正準形）でない状態式で与えられる場合には、それを上記 へ変換しないと(7-19)のような方法が適用できない。この場合には次のAckermann公式が役立つ 任意の行列Pで代表されるｎ次プラント 𝒙 𝑘 + 1 = 𝑷𝒙 𝑘 + 𝒒𝒖 𝑘 において、状態ベクトルFeedback, u(k)=-kx(k) のゲインｋは次のように決定する 𝒌= 0 0 ⋯ 1 𝒒 𝑷𝒒 ⋯ 𝑷𝑛−1 𝒒 −1 𝐹(𝑷) (7-20) ただし、係数 α1, α1,…, αn が式(7-15)のように一組の指定固有値を代表する 𝐹 𝑷 = 𝑷𝑛 + 𝛼1 𝑷𝑛−1 + ⋯ + 𝛼𝑛 𝑰 (7-21) I 動作つき状態ベクトルフィードバック系の場合にはPの中にI 動作も含ませる。したがって、 (n+1)次系となり、(7-21)の係数が(7-15)とおなじくαn+1 に達する

7.3 固有値を指定の状態推定系設計 プラント: 𝒙 𝑘 + 1 = 𝑷𝒙 𝑘 + 𝒒𝒖 𝑘 𝒚 𝑘 = 𝒄𝒙(𝑘) についての状態推定アルゴリズム（観測器）は 28 ここからは状態推 定系について 暫定推定値 step1: ෝ° 𝑘 + 1 = 𝑷ෝ 𝒙 𝒙 𝑘 + 𝒒𝑢(𝑘) step2: ෝ 𝑘+1 =𝒙 ෝ° 𝑘 + 1 + 𝒇 𝑦 𝑘 + 1 − 𝒄ෝ 𝒙 𝒙° 𝑘 + 1 両式をまとめると ෝ 0 =𝟎 𝒙 (7-22) ෝ 𝑘 + 1 = 𝑷ෝ 𝒙 𝒙 𝑘 + 𝒒𝑢(𝑘) + 𝒇𝒄 𝒙 𝑘 + 1 − 𝑷ෝ 𝒙 𝑘 − 𝒒𝑢(𝑘) ෝ 𝑘 となるから推定誤差ベクトル 𝒆 𝑘 = 𝒙 𝒌 − 𝒙 は、 𝒆 𝑘 + 1 = 𝑷 − 𝒇𝒄𝑷 𝒆 𝑘 (7-23) となり、制御系から完全に分離した動的システムを形成する。しかも 𝒙 𝑘 + 1 = 𝑷 − 𝒒𝒌 𝒙 𝑘 により表せれる状態ベクトルフィードバック制御系との間に、次の対応関係をもつ

7.3 固有値を指定の状態推定系設計 29 状態ベクトルフィードバック制御系と状態推定系の間に、次の対応関係がある 制御系 推定系 𝑷 − 𝒒𝒌 𝑷 − 𝒇𝒄𝑷 𝒌 𝒇T 𝒒 𝒄𝑷 𝑷 𝑷T T この双対関係に注目すると式(7-20)に対応する公式は 𝒇 = 𝐹(𝑷) ただし 𝒄 𝒄𝑷 𝒄𝑷𝟐 ⋮ 𝒄𝑷𝑛−1 −1 0 0 ⋮ 0 1 (7-24) 𝐹 𝑷 = 𝑷𝑛−1 + 𝛼1 𝑷𝑛−2 + ⋯ + 𝛼𝑛−1 𝑷 + 𝛼𝑛 𝑰 として (n-1) の固有値を指定し、第nの固有値は０とする。式(7-24)を使うと式 (7-23)が代表する観測系へ (n-1) 個の固有値を付与するようなfが決定される

7.4 有限整定系 7.1 LQ制御系の固有値 7.2 LQ制御系の根軌跡 7.3 固有値を指定の制御系および状態推定系設計 7.4 有限整定系 7.5 多変数制御系 すべての固有値が0のプラントの応答は有限整定する ⚫ プロセス特性Gp(z)に時系列形(実験 データ gi)を用いたとき、その giが 少しぐらい違ってもに深刻な影響は ない 30

7.4 有限整定系とは 31 有限時間整定制御（デッドビート制御） 通常の線形制御では出力応答は指数関数的な挙動となり時刻∞で目標値に達するような応答になりま す。有限時間整定制御は，有限の指定された時刻で目標値に達し，以後，目標値と一致し続けるよう な制御です。離散時間系の有限整定，連続時間系の有限整定（むだ時間を使用する），非線形同次シ ステムを用いた有限整定などがあります。また，終端状態制御も有限整定制御の一種と捉えることが できます。 https://www.sice-ctrl.jp/about/ctrl-info/detail/ 例）ディジタル制御 その4 v2 6.5 ０型プラントのLQI制御(observer出力) 推定値(青)とその状態値(赤)を重ねて表示 n=３のシステム 有限整定観測では３刻みで収束している （observer初期値は異なる値を設定した） 対象システム： 𝐺 𝑧 = DGC_no4/p118v1 B 𝑧 𝑏0 𝑧 + 𝑏1 = 2 𝐴 𝑧 𝑧 + 𝑎1 𝑧 + 𝑎2

• https://www.sice-ctrl.jp/about/ctrl-info/detail/

7.4 有限整定系(例1） 例として純むだ時間プラントを取り上げる： 例）n=3 step1：可制御正準形としてP, q, cを求める 1 𝐺𝑝 = 3 𝑧 可制御 正準形 𝑏1 = 0, 𝑏2 = 0, 𝑏3 = 1 𝑎1 = 0, 𝑎2 = 0, 𝑎3 = 0 0 𝑃= 0 −𝑎3 1 0 −𝑎2 𝑦 𝑘+1 𝒔 𝑘+1 = 0 𝑞= 0 1 0 1 −𝑎1 1 𝒄𝑷 𝑦 𝑘 0 𝑷 𝒔 𝑘 𝑷1 ∈ ℝ 𝑐 = 𝑏3 ０型プラントのLQI制御 𝒒1 ∈ ℝ 𝑏2 𝑏1 (3-14) 𝒙 𝑘 + 1 = 𝑷𝒙 𝑘 + 𝒒𝑢(𝑘) 𝑦 𝑘 = 𝒄𝒙(𝑘) 𝒄𝒒 + 𝒒 𝑑(𝑘) 𝑛+1 × 𝑛+1 step３: LQI制御プログラム実行 既出：6.5 1 𝐺𝑝 = 𝑛 𝑧 B(z) 𝑏1 𝑧 𝑛−1 + 𝑏2 𝑧 𝑛−2 + ⋯ + 𝑏𝑛−1 𝑧 + 𝑏𝑛 𝐺 𝑧 = = A(z) 𝑧 𝑛 + 𝑎1 𝑧 𝑛−1 + 𝑎2 𝑧 𝑛−2 + ⋯ + 𝑎𝑛−1 𝑧 + 𝑎𝑛 step2： LQI制御系へ変換（P1,q1を算定する） 新たな 状態方程式： 32 𝑛+1 ×1 (6-24)

7.4 有限整定系（例1） 33 １．有限整定系におけるLQI制御での状態ベクトルフィードバックGain 1 𝐺𝑝 = 𝑛 ： これの状態ベクトルフィードバックGainは次のようになる 𝑧 𝐾0 = 𝐾2 = 𝐾2= ⋯= 𝐾𝑛 = 𝐾, 𝐾1 = 0 𝐾= (7-26) 2 1 + 1 + 4𝑤 ２．次に根についてみてみる I動作を直列結合した開ループ伝達関数： LQI制御系の特性式： 𝐹𝑐 𝑧 = 𝐴 𝑧 ∙ 𝐴 𝑧 −1 + 𝐺 𝑧 = 𝑧 1 𝐵(𝑧) = 𝑧 − 1 𝑧 𝑛 𝐴(𝑧) 1 𝐵 𝑧 ∙ 𝐵 𝑧 −1 𝑤 1 1 𝐹𝑐 𝑧 = 𝑧 − 1 𝑧 𝑛 𝑧 −1 − 1 𝑧 −𝑛 + 𝑤 𝑧𝑧 −1 = 𝑧 − 1 𝑧 −1 − 1 + 𝑤 1 = 𝑧 −1 − 𝑧 + 2 + 𝑤 =0 根は次式の解として求まる 𝑧2 + 2 + 1 𝑧+1=0 𝑤

7.4 有限整定系（ 例1 Step応答） 純むだ時間プラントのStep応答 34 純むだ時間系のLQI制御例 𝐺𝑝 = 1 𝑧3 step状外乱 (step:25で) V=0.5を入力 伝達関数の通りに3刻み遅延 Gain= [[0.9161, 0.0 , 0.9161, 0.9161]] 𝐾= noise 2 1 + 1 + 4𝑤 = 2 1 + 1 + 4 ∗ 0.1 = 0.91608 根= [11.9161 0.0839] 不安定根 𝑧2 1 + 2+ 𝑧+1=0 0.1 DGC_no5/ p134_n3

7.4 有限整定系 35 以上のような簡単さは、制御系の設計段階でも実施についても、いろいろな示唆 を与える。たとえば設計のさいに、まず可制御標準形のプラントモデル 𝑘1 = −𝑎𝑛 , 𝑘2 = −𝑎𝑛−1 , ⋯, 𝑘𝑛 = −𝑎1 の状態ベクトルフィードバックを想定して、プラントを式(7-25)のむだ時間系 に変える。次いで式(7-26)のKを適当にえらんだ制御をつけ加える。これはもは や最適制御ではなく、式(6-20)のような評価関数の含む項wd2(k)は好ましい制 御を得るためのめやすにすぎない。したがって設計者は評価関数値の最適化にこ だわりすぎる必要はない。操作部にヒステリシスのような非線形特性があれば、 どのみち最適制御は通用しなくなる。 （原文そのまま） “高橋安人,ディジタル制御 ,岩波書店,1987, p135 1 𝐺𝑝 = 𝑛 𝑧 ∞ 𝐾0 = 𝐾2 = 𝐾2 = ⋯= 𝐾𝑛 = 𝐾, 𝐾1 = 0 (7-25) 𝐾= 𝐽 = ෍ 𝑦 2 𝑘 + 1 + 𝑤𝑑 2 𝑘 2 1 + 1 + 4𝑤 𝑘=0 (7-26) (6-20)

7.4 有限整定系（例２：時系列データ） 例２：（時系列データ） プロセス特性Gp(z)に時系列形(実験データ gi)を 用いて有限制御Gainを求めたとき、現実のプラ ント（固定）と推定パラメータ（変化させる）と の不一致の影響を確かめる ±１０％、２０％、 ３０％変化 例で利用するプラントモデル G 𝑧 = 𝑧 𝑏1 𝑧 + 𝑏2 𝐵(𝑧) = 𝑧 − 1 𝑧 − 𝑝1 𝑧 − 𝑝2 𝐴(𝑧) （4.3節参照） 36

7.4 有限整定系（例２：時系列データ） 37 １．プロセス特性Gp(z)に時系列形(実験データ gi)を用いる 𝐺𝑝 = 𝑔1 𝑧 −1 + 𝑔2 𝑧 −2 + 𝑔3 𝑧 −3 ⋯ + 𝑔𝑛 𝑧 −𝑛 Τ 1 − 𝑝𝑧 −1 𝑔1 , 𝑔2 , ⋯ , 𝑔𝑛 実測した応答データ ２．状態方程式を求める 0 0 0 𝑃= ⋮ 0 0 0 伝達関数形 1 0 0 ⋮ 0 0 0 0 0 1 0 0 1 ⋯ ⋯ ⋯ 0 ⋯ ⋯ ⋯ ⋯ 𝒙 𝑘 + 1 = 𝑷𝒙 𝑘 + 𝒒𝒖 𝑘 ⋯ ⋯ ⋯ ⋯ 1 0 0 0 0 0 ⋮ 0 𝑔6 𝑝 𝑔1 𝑔2 𝑔3 𝑞= ⋮ ⋮ 𝑔𝑛−1 1 𝐶= 1 𝑏1 𝑧 𝑛−1 + 𝑏2 𝑧 𝑛−2 + ⋯ + 𝑏𝑛−1 𝑧 + 𝑏𝑛 𝐵(𝑧) 𝐺𝑝 = = 𝑧 𝑛−1 𝑧 − 𝑝 𝐴(𝑧) （4.3節参照） 𝒚 𝑘 = 𝒄𝒙(𝑘) 0 ⋯ ⋯ 0 0 𝑏1 = 𝑔1 , 𝑏2 = 𝑔2 − 𝑝𝑔1 , ⋯ , 𝑏𝑛 = 𝑔𝑛 − 𝑝𝑔𝑛−1 (7-27) (7-28) ３．この有限整定制御Gainは次のようになる 𝐾0 = 1Τ𝐵 1 = 1Τ 𝑏1 + 𝑏2 + ⋯ + 𝑏𝑛 , 𝐾1 = 0 𝐾2 = 𝐾2 = ⋯= 𝐾𝑛−1 = 𝐾0 𝐾𝑛 = 1 + 𝑝 − 𝐾0 𝑔1 + ⋯ + 𝑔𝑛−1 (7-29)

7.4 有限整定系(例２：時系列データ） 38 １．プロセス特性Gp(z)に時系列形(実験データ gi)を用いる ここでは、本末転倒かもしれないが計算で求める 対象プロセス： 𝐾𝑒 −𝑠𝐿 𝐺 𝑠 = 1 + 𝑇1 𝑠 1 + 𝑇2 𝑠 L：遅延時間 Step1: G(s)→G(z)変換する（方法１ or 方法２） Step2: G(z)を用い計算によって、g1,g2,g3…を求める （4.3節参照） 計算で時系列形 (実験データ gi)を 求める G(s)→G(z)変換方法として2つある 方法１：状態空間における数値計算 （４．２節） 方法２：解析的方法 (4-12) 式 （ ４．３節） 今回、方法２を用いてG(s)→G(z)変換を行ったところ、書籍表7-1と若干異なる 結果となった。そこで、方法１で数値計算により求めてみた。すると結果は、方 法１と方法２で算出した値は等しくなった。これより、書籍表7-1の値は、もしか したら違っている可能性がある。よって、表7-1に関わる計算は、新たに計算した 値を利用することにする。

7.4 有限整定系（表7-1） 書籍の表7-1と方法2での数値は、若干異なる 値となった そこで、方法１で数値計算により求めてみた。す ると、方法１と方法２での計算結果は一致した。 これより、書籍表7-1の値は、もしかしたら違って いる可能性がある。よって、表7-1に関わるシミュ レーションには、新たに計算した値を利用するこ （私の計算が間違っている可能性も高い） とにした。 39

40 方法１：状態空間における数値計算（４.２節） case (b) L1≠0 𝐺 𝑧 = 𝐵1 𝑧 𝐵2 𝑧 −1 + 𝑧 𝐴 𝑧 𝐴 𝑧 𝑄1 𝑄2 u(k) u(k-1) L1 (k-1)T 𝑸0 = න 𝑒 𝑨𝜂 𝑑𝜂 𝑄2 𝑇−𝐿1 𝑸1 = න 0 (k+1)T 伝達関数はQ1とQ2 の和になる t 𝑸1 𝑇 T-L1 kT 1個時間遅れ入力 𝑸0 連続時間プラントの入出力間に むだ時間遅れL1があると、階段 状入力が右図のようにズレる。 y(k+1) y(k) 𝑄1 A(z)：全て同じ 𝑸2 𝑒 𝑨𝜂 𝑑𝜂 𝑇 𝑒 𝑨𝜂 𝑑𝜂 = 𝑸0 − 𝑸1 𝑸2 = න 0 𝑇−𝐿1 (list 4-1) p57v1.py 𝒒 = 𝑸0 𝒃 (list 3-1) 𝒒𝟏 = 𝑸1 𝒃 𝒒𝟐 = 𝑸2 𝒃 page40v1.py 𝑷𝒄 = 𝑻𝑷𝑻−1 𝒒𝒄 = 𝑻𝒒 𝑪𝒄 = 𝑪𝑻−1 𝑷𝒄 = 𝑻𝑷𝑻−1 𝒒𝒄𝟏 = 𝑻𝒒1 𝑪𝒄𝟏 = 𝑪𝑻−1 𝑷𝒄 = 𝑻𝑷𝑻−1 𝒒𝒄𝟐 = 𝑻𝒒2 𝑪𝒄𝟐 = 𝑪𝑻−1 𝐴 = 𝑎𝑛 𝑎𝑛−1 𝑎𝑛−2 ⋯ 𝑎 𝐴 = 𝑎𝑛 𝑎𝑛−1 𝑎𝑛−2 ⋯ 𝑎 𝐴 = 𝑎𝑛 𝐵 = 𝑏𝑛 𝑏𝑛−1 𝑏𝑛−2 ⋯ 𝑏1 𝐵1 = 𝑏𝑛 𝑏𝑛−1 𝑏𝑛−2 ⋯ 𝑏1 𝐵2 = 𝑏𝑛 u(k) 𝑎𝑛−1 𝑎𝑛−2 ⋯ 𝑏𝑛−1 𝑏𝑛−2 ⋯ 𝑏1 u(k-1) 𝑎

41 方法２：解析的方法（ ４.３節） むだ時間プロセス系のG(s)をG(z)に変換した例を示す （プラントの説明は（A-3)を参照のこと） パラメータ変換表 𝐾𝑒 −𝑠𝐿 𝐺 𝑠 = 1 + 𝑇1 𝑠 1 + 𝑇2 𝑠 (4-10) (4-12) 𝑇1 ≠ 𝑇2 部分分数展開→表4-4 を利用→まとめる 𝑘 𝑏1 𝑧 2 + 𝑏2 𝑧 + 𝑏3 𝐺 𝑧 = 𝑁+1 2 𝑧 𝑧 + 𝑎1 𝑧 + 𝑎2 遅れ時間Lとサンプリング時間Tの関係 (4-11) L=NT+L1 （遅れ時間がTより大きくなるとN>=１） 圧力、拡散、混合、反応過程に、しばしば(4-10)ような例はあるそうだ。このため工業プロセスではステップ応答の 実験結果に合うように(4-10)式のK,L,T1,T2を決めることがあるとのこと。流体力学、拡散、反応などを含むプロセ スは複雑であり、モデルをシステム解析から求めるより、実験的に(4-10)のパラメータを見出す方が近路とのことが 多い。（原著p237）

42 シミュレーションに使う時系列形データ gi gi (a) (b) 表7-1 G(z)のパラメータ 真のプラントよりT1,T2,L1を±10％変化させる ① (c) a,bパラメータ case(a)：書籍と計算結果一致 case(b) ：書籍と計算わずか異なる case(c) ：書籍と計算わずか異なる ② ｐパラメータは全て一致した #case(a) i= 6, 0.70788 #case(b) i= 6, 0.73687 #case(c) i= 5, 0.68864 DGC_no5/ p136a

7.4 有限整定系（表7-1） 43 ｐ.136-137 G(s) G(z) “高橋安人,ディジタル制御 ,岩波書店,1987

T1,T2,L1を２０％変化 (d) ２０％大きくする (e) ２０％小さくする 44 (f) ３０％大きくする T1,T2,L1を３０％変化 シミュレーションに使う時系列形データ gi (g) ３０％小さくする

7.4 有限整定系(例２：時系列データ） 45 T1,T2,Lは変化なし 実際のプロセス（固定） p0=0.708;g1=0.008;g2=0.128;g3=0.183; g4=0.172;g5=0.141;g6=0.108 #case(a) 制御に用いた値 p0=0.708;g1=0.008;g2=0.128;g3=0.183; g4=0.172;g5=0.141;g6=0.108 Gain= [[ 3.41829 0.00000 3.41829 3.41829 3.41829 3.41829 0.45236]] noise F= [[1.00000] [0.70800] [0.50126] [0.35489] [0.25127] [0.17790] [0.12595]] step状外乱 (step:25で) V=0.5を入力 yopen： 開ループステップ応答 (状態フィードバックなし) DGC_no5/ p138

7.4 有限整定系(例２ ：時系列データ） 46 10%大きい 実際のプロセス（固定） p0=0.708;g1=0.008;g2=0.128;g3=0.183; g4=0.172;g5=0.141;g6=0.108 #case(a) 制御に用いた値 p0=0.737;g1=0.004;g2=0.112;g3=0.179; g4=0.179;g5=0.155;g6=0.124 #case(b) Gain= [[ 3.80046 0.00000 3.80046 3.80046 3.80046 3.80046 -0.43578]] noise F= [[1.00000] [0.73687] [0.54298] [0.40010] [0.29482] [0.21725] [0.16008]] step状外乱 (step:25で) V=0.5を入力 DGC_no5/ p138

7.4 有限整定系(例２：時系列データ） 47 10%小さい 実際のプロセス（固定） p0=0.708;g1=0.008;g2=0.128;g3=0.183; g4=0.172;g5=0.141;g6=0.108 #case(a) 制御に用いた値 p0=0.6886;g1=0.01531;g2=0.15946;g3= 0.20368;g4=0.17896;g5=0.1378#;g6=0. 09935 #case(c) Gain= [[ 3.03777 0.00000 3.03777 3.03777 3.03777 3.03777 -0.42329]] noise F= [[1.00000] [0.68860] [0.47417] [0.32651] [0.22484] [0.15482] [0.10661]] step状外乱 (step:25で) V=0.5を入力 DGC_no5/ p138

７.4 有限整定系(例２：時系列データ） 48 20%大きい 破綻？ 実際のプロセス（固定） p0=0.708;g1=0.008;g2=0.128;g3=0.183; g4=0.172;g5=0.141;g6=0.108 #case(a) 制御に用いた値 p0=0.7619;g1=0.0009;g2= 0.08062;g3= 0.14412;g4= 0.15281;g5= 0.13792;g6= 0.11514 #case(f1) noise Gain= [[ 4.20013 0.00000 4.20013 4.20013 4.20013 4.20013 0.40692]] F= [[1.00000] [0.76190] [0.58049] [0.44228] [0.33697] [0.25674] [0.19561]] DGC_no5/ p138

７.4 有限整定系(例２：時系列データ） 49 20%小さい 実際のプロセス（固定） p0=0.708;g1=0.008;g2=0.128;g3=0.183; g4=0.172;g5=0.141;g6=0.108 #case(a) 制御に用いた値 p0=0.6477;g1=0.02797;g2=0.19767;g3=0.2 2438;g4=0.1815;g5=0.12982;g6=0.08732 #case(g1) noise Gain= [[ 2.81262 0.00000 2.81262 2.81262 2.81262 2.81262 0.49366]] F= [[1.00000] [0.64770] [0.41952] [0.27172] [0.17599] [0.11399] [0.07383]] DGC_no5/ p138

7.4 有限整定系(例２：時系列データ） 50 30%大きいい 破綻 実際のプロセス（固定） p0=0.708;g1=0.008;g2=0.128;g3=0.183; g4=0.172;g5=0.141;g6=0.108 #case(a) 制御に用いた値 p0=0.7723;g1=0.00005;g2=0.06294;g3= 0.12728;g4=0.14207;g5=0.13349;g6=0. 11552;g7=0.09534#case(d1) noise Gain= [[ 4.51285 0.00000 4.51285 4.51285 4.51285 4.51285 -0.32992]] F= [[1.00000] [0.77230] [0.59645] [0.46064] [0.35575] [0.27475] [0.21219]] DGC_no5/ p138

7.4 有限整定系(例２：時系列データ） 51 30%小さい 実際のプロセス（固定） p0=0.708;g1=0.008;g2=0.128;g3=0.183; g4=0.172;g5=0.141;g6=0.108 #case(a) 制御に用いた値 p0=0.5987;g1= 0.04844; g2=0.2435; g3=0.24235; g4=0.17731; g5=0.11573; g6=0.07138 #case(g1) noise Gain= [[ 2.47901 0.00000 2.47901 2.47901 2.47901 2.47901 0.45226]] F= [[1.00000] [0.59870] [0.35844] [0.21460] [0.12848] [0.07692] [0.04605]] step状外乱 (step:25で) V=0.5を入力 DGC_no5/ p138

7.4 有限整定系(パラメータ不一致の影響） 52 ⚫ プラントと設計に用いたパラメータが一致しな いために、これらの計算結果は真の有限整定で はない。 ⚫ しかし、この例に関する限りは、１０％の不一 致はさして深刻ではない。 ⚫ ただし、２０～３０％程度違ったときには、な んらかの不具合をみることができる。

7.4 有限整定系（例２ ：時系列データ） “高橋安人,ディジタル制御 ,岩波書店,1987 53 書籍の図 page139 ? (7-29)式と異なる 𝐾𝑛−1 = 1 + 𝑝 − 𝐾0 𝑔1 + ⋯ + 𝑔𝑛−1 (7-29)式による計算結果： [[ 3.45510 0.00000 3.45510 3.45510 3.45510 3.45510 -0.43626]] (7-29)

7.4 有限整定系（例２：observer） 54 １．プロセス特性Gp(z)に時系列形(実験データ gi)を用いる 𝐺𝑝 = 𝑔1 𝑧 −1 + 𝑔2 𝑧 −2 + 𝑔3 𝑧 −3 ⋯ + 𝑔𝑛 𝑧 −𝑛 Τ 1 − 𝑝𝑧 −1 𝑔1 , 𝑔2 , ⋯ , 𝑔𝑛 実測した応答データ ２．状態方程式を求める 0 0 0 𝑃= ⋮ 0 0 0 伝達関数形 1 0 0 ⋮ 0 0 0 0 0 1 0 0 1 ⋯ ⋯ ⋯ 0 ⋯ ⋯ ⋯ ⋯ 𝒙 𝑘 + 1 = 𝑷𝒙 𝑘 + 𝒒𝒖 𝑘 ⋯ ⋯ ⋯ ⋯ 1 0 0 0 0 0 ⋮ 0 𝑔6 𝑝 𝑔1 𝑔2 𝑔3 𝑞= ⋮ ⋮ 𝑔𝑛−1 1 𝐶= 1 （4.3節参照） 𝒚 𝑘 = 𝒄𝒙(𝑘) 0 ⋯ 𝑏1 𝑧 𝑛−1 + 𝑏2 𝑧 𝑛−2 + ⋯ + 𝑏𝑛−1 𝑧 + 𝑏𝑛 𝐵(𝑧) 𝐺𝑝 = = 𝑧 𝑛−1 𝑧 − 𝑝 𝐴(𝑧) ⋯ 0 (7-27) 0 𝑏1 = 𝑔1 , 𝑏2 = 𝑔2 − 𝑝𝑔1 , ⋯ , 𝑏𝑛 = 𝑔𝑛 − 𝑝𝑔𝑛−1 (7-28) ３．Observerに用いる出力フィードバックゲイン行列fは次のようになる 𝑓𝑇 = 1 𝑝 𝐹 𝑇 = 𝑰 𝑹 𝑹2 𝑝2 ⋯ 𝑝𝑛−1 ⋯ 𝑹𝑛−2 𝑹𝑛−1 𝑮−1 𝑛 （7.4 次頁） （7.5 多変数系で利用） (7-30)

例２：observer出力 T1,T2,Lは変化なし 実際のプロセス（固定） p0=0.708;g1=0.008;g2=0.128;g3=0.183; g4=0.172;g5=0.141;g6=0.108 #case(a) 制御に用いた値 p0=0.708;g1=0.008;g2=0.128;g3=0.183; g4=0.172;g5=0.141;g6=0.108 Gain= [[ 3.41829 0.00000 3.41829 3.41829 3.41829 3.41829 0.45236]] F= [[1.00000] [0.70800] [0.50126] [0.35489] [0.25127] [0.17790] [0.12595]] step状外乱 (step:25で) V=0.5を入力 55 DGC_no5/ p138

7.5 多変数制御系 7.1 LQ制御系の固有値 7.2 LQ制御系の根軌跡 7.3 固有値を指定の制御系および状態推定系設計 →1入力1出力（SISO） 7.4 有限整定系 7.5 多変数制御系 →2入力２出力（MIMO） ⚫ 7.4の多変数制御系への拡張 ⚫ 多変数間の干渉の問題 56

7.5 多変数制御系 ｍ入力ｍ出力（MIMO）（2≤m） 1入力1出力（SISO） 0 0 0 𝑃= ⋮ 0 0 0 1 0 0 ⋮ 0 0 0 0 0 1 0 0 1 ⋯ ⋯ ⋯ 0 ⋯ ⋯ ⋯ ⋯ ⋯ 0 ⋯ 0 ⋯ 0 ⋯ ⋮ 1 0 0 𝑔6 0 𝑝 𝑔1 𝑔2 𝑔3 𝑞= ⋮ ⋮ 𝑔𝑛−1 1 𝐶= 1 0 ⋯ ⋯ 0 0 (7-27) 𝐺𝑝 = 𝑔1 𝑧 −1 + 𝑔2 𝑧 −2 + ⋯ + 𝑔𝑛 𝑧 −𝑛 Τ 1 − 𝑝𝑧 −1 𝑏1 = 𝑔1 , 𝑏2 = 𝑔2 − 𝑝𝑔1 , ⋯ , 𝑏𝑛 = 𝑔𝑛 − 𝑝𝑔𝑛−1 (7-28) 𝐾𝑛 = 1 + 𝑝 − 𝐾0 𝑔1 + ⋯ + 𝑔𝑛−1 𝟎 𝟎 𝟎 𝑷= ⋮ 𝟎 𝟎 𝟎 𝑰 𝟎 𝟎 ⋮ 𝟎 𝟎 𝟎 𝑪= 𝑰 𝟎 𝟎 𝟎 𝑰 𝟎 𝟎 𝑰 ⋯ ⋯ ⋯ 𝟎 ⋯ ⋯ ⋯ ⋯ (7-26) ⋯ ⋯ ⋯ ⋯ 𝑰 𝟎 𝟎 ⋯ ⋯ 𝟎 𝟎 𝟎 𝟎 ⋮ 𝟎 𝑮𝒏 𝑹 𝑮𝟏 𝑮𝟐 𝑮𝟑 Q= ⋮ ⋮ 𝑮𝒏−𝟏 𝑰 𝟎 𝑮 𝑧 = 𝑩𝟏 𝑧 𝑛−1 + ⋯ + 𝑩𝒏 ∙ 𝑧 𝑛−1 𝑧𝑰 − 𝑹 (4-22) −1 𝑩𝟏 = 𝑮𝟏 , 𝑩𝟐 = 𝑮𝟐 − 𝑮𝟏 𝑹, ⋯ , 𝑩𝒏 = 𝑮𝒏 − 𝑮𝒏−𝟏 𝑹 (m=2のとき） 𝑝 𝑹= 1 0 𝑔 (𝑖) 𝑔12 (𝑖) 0 , 𝑮𝒊 = 11 𝑝2 𝑔21 (𝑖) 𝑔22 (𝑖) 𝑲𝟎 = 𝑩𝟏 + ⋯ + 𝑩𝒏 𝐾0 = 1Τ𝐵 1 = 1Τ 𝑏1 + 𝑏2 + ⋯ + 𝑏𝑛 , 𝐾1 = 0 𝐾2 = 𝐾2 = ⋯= 𝐾𝑛−1 = 𝐾0 57 −𝟏 , (4-21) 𝑲𝟏 = 𝟎 𝑲𝟐 = 𝑲𝟐 = ⋯= 𝑲𝒏−𝟏 = 𝑲𝟎 𝑲𝒏 = 𝟏 + 𝑹 − 𝑲𝟎 𝑮𝟏 + ⋯ + 𝑮𝒏−𝟏 (7-31)

7.5 多変数制御系 58 1. １入出力数の有限整定制御を多入出力系へ拡張する 𝑮 𝑧 = 𝑩𝟏 𝑧 𝑛−1 + ⋯ + 𝑩𝒏 ∙ 𝑧 𝑛−1 𝑧𝑰 − 𝑹 −1 (4-21) 𝑩𝟏 = 𝑮𝟏 , 𝑩𝟐 = 𝑮𝟐 − 𝑮𝟏 𝑹, ⋯ , 𝑩𝒏 = 𝑮𝒏 − 𝑮𝒏−𝟏 𝑹 例）m=2,n=5 実験の 時系列 データ 𝑹= 𝑝1 0 𝑔 (𝑖) 0 , 𝑮𝒊 = 11 𝑝2 𝑔21 (𝑖) 𝑔12 (𝑖) , i = 1,2, ⋯ , 𝑛 𝑔22 (𝑖) 𝑢1 → 𝑦1 : 𝑔11 (0), 𝑔11 (1), 𝑔11 (1) , 𝑔11 (3) , 𝑔11 (4) 𝑢1 → 𝑦2 : 𝑔12 (0),𝑔12 (1), 𝑔12 (1) , 𝑔12 (3) , 𝑔12 (4) 𝑢2 → 𝑦1 : 𝑔21 (0), 𝑔21 (1), 𝑔21 (1) , 𝑔21 (3) , 𝑔21 (4) 𝑢2 → 𝑦2 : 𝑔22 (0),𝑔22 (1), 𝑔22 (1) , 𝑔22 (3) , 𝑔22 (4) Rが対角形であるのは応答に含まれる等比を入力だけから決めることを意味する。たとえばu1→y1, u1→y2の 応答成分がともに等比p1をもって減衰するものと考える。 𝟎 𝟎 𝟎 𝑷= ⋮ 𝟎 𝟎 𝟎 𝑰 𝟎 𝟎 ⋮ 𝟎 𝟎 𝟎 𝟎 𝟎 𝑰 𝟎 𝟎 𝑰 ⋯ ⋯ ⋯ 𝟎 ⋯ ⋯ ⋯ ⋯ ⋯ ⋯ ⋯ ⋯ 𝑰 𝟎 𝟎 𝟎 𝟎 𝟎 ⋮ 𝟎 𝑮𝒏 𝑹 𝑮𝟏 𝑮𝟐 𝑮𝟑 Q= ⋮ ⋮ 𝑮𝒏−𝟏 𝑰 𝑪= 𝑰 𝟎 ⋯ ⋯ 𝟎 𝟎 (4-22) 𝑷 ∈ ℝ𝑚𝑛×𝑚𝑛 , 𝑸 ∈ ℝ𝒎𝑛×𝑚 , 𝑪 ∈ ℝ𝑚×𝑚𝑛

7.5 多変数制御系 59 2. １入出力数の有限整定制御を多入出力系へ拡張する このm変数プロセスの有限整定行列は式(7-19)と同じ形をした次式になる 𝑲𝟎 = 𝑩𝟏 + ⋯ + 𝑩𝒏 −𝟏 𝑲𝟏 = 𝟎 (7-31) 𝑲𝟐 = 𝑲𝟐 = ⋯= 𝑲𝒏−𝟏 = 𝑲𝟎 𝑲𝒏 = 𝟏 + 𝑹 − 𝑲𝟎 𝑮𝟏 + ⋯ + 𝑮𝒏−𝟏 dc-gain: 𝑮 𝟏 = 𝑪 𝑰−𝑷 −𝟏 𝑸 III=np.eye(10,10) invP=np.linalg.inv(III-P) print(np.dot(C,np.dot(invP,Q)))

7.5 多変数制御系（例：2変数プロセス） １．プロセス特性Gp(z)に時系列形(実験データ gi)を用いる 例）2変数プロセス For input u1: p1=0.8 2変数プロセスモデルを作成 For input u2: p2=0.6 u1→y1 g(1)=0.1 g(2)=0.21 g(3)=0.18 g(4)=0.15 g(5)=0.13 u2→y1 g(1)=0.02 g(2)=0.06 g(3)=0.14 g(4)=0.18 g(5)=0.12 u1→y2 g(1)=-0.05 g(2)=-0.04 g(3)=0.01 g(4)=0.1 g(5)=0.09 u2→y2 g(1)=0.08 g(2)=0.2 g(3)=0.23 g(4)=0.19 g(5)=0.14 時系列形(実験データ gi) 𝑹= 𝑔 (𝑖) 0 , 𝑮𝒊 = 11 𝑝2 𝑔21 (𝑖) 𝑝1 0 𝟎 𝟎 𝟎 𝑷= ⋮ 𝟎 𝟎 𝟎 𝑰 𝟎 𝟎 ⋮ 𝟎 𝟎 𝟎 𝟎 𝑰 𝟎 ⋯ ⋯ ⋯ ⋯ 𝟎 𝟎 𝑰 ⋯ 𝟎 ⋯ ⋯ 𝑪= 𝑰 𝟎 ⋯ ⋯ ⋯ ⋯ 𝑰 𝟎 𝟎 ⋯ 𝟎 𝟎 𝟎 ⋮ 𝟎 𝑮𝒏 𝑹 ⋯ 𝑔12 (𝑖) 𝑔22 (𝑖) 𝑮𝟏 𝑮𝟐 𝑮𝟑 Q= ⋮ ⋮ 𝑮𝒏−𝟏 𝑰 𝟎 𝟎 60

7.5 多変数制御系（例：2変数プロセス） 61 ステップ入力に対する応答 ⚫ u1→ y2、 u2→ y1という 入出力関係がかなり強い u1--->y1, y2 u2--->y1, y2 ⚫ このため、図7-7の刻み10 で外乱v1が加わるとy1ばか りでなく、y2も乱される ⚫ 同様に図7-7の刻み２0で 外乱v2によりy1も乱れる dc-gainを計算してみる 1. dc-gainの計算結果、 すなわち時間無限大で の応答を示す 2. シミュレーションstep 応答の最終値と合って いる In=np.eye(n,n) 𝑮 𝟏 = 𝑪 𝑰−𝑷 −𝟏 𝑸 3. 但し、書籍の値と違う invIP=np.linalg.inv((In-P)) dcGain=np.dot(C,np.dot(invIP,Q)) print("dc_Gain=¥n",dcGain) DGC_no5/p142DCG

7.5 多変数制御系（例：2変数プロセス） 時系列形(実験データ gi) 62 p141 p141のdc-gain値: 図7-6 step応答最終値と合っている （注：前頁と表示順序と違う） ? p142 “高橋安人,ディジタル制御 ,岩波書店,1987

7.5 多変数制御系（例：2変数プロセス） 63 有限整定制御 干渉 干渉 noise noise step外乱 (K=20でv1, k=20でv2) V=1を入力 DGC_no5/ p142v1

7.5 多変数制御系（例：2変数プロセス） 64 p142 “高橋安人,ディジタル制御 ,岩波書店,1987

7.5 多変数制御系（例：2変数プロセス） 65 1. このシミュレーションでは10刻み毎に設定点r1, r2 をステップ状に変化させた り、ステップ状外乱v1, v2を加えた。 2. 図７－６に示したようにu1--->y2, u2--->y1という入出力関係がかなり強いの で、図7-7の刻み１０で外乱v1が加わるとy1ばかりでなくy2も乱される。刻み ２０でも外乱v2が加わるとy2も乱される。 3. 多変数制御で大きな話題とされてきたのは、設定点入力についての非干渉化であ る。 4. 刻み３０でy2の設定点を１--->2へ上げるとy1がふれる。また刻み40では、y1 の設定点が３--->４と上げるとy2の値を乱している。 5. このような、干渉をできるだけ除くような方法として、主ループへゲイン行列 (statics decoupler)を挿入する方法が知られている。非干渉化はプラントの行 列伝達関数の対角化である。それを完全な対角化ではなく、対角項が優位 (diagonal dominance) になるような近似がある。次頁で記述する （原文をすこし省略改変）

7.5 多変数制御系（非干渉化） 非干渉化：行列伝達関数の対角項が優位になるように近似する ∞ 評価関数： 𝐽=෍ 𝒚 𝑘 + 1 − 𝒚𝑓 𝑇 𝑾𝑦 𝒚 𝑘 + 1 − 𝒚𝑓 + 𝒖 𝑘 − 𝒖𝑓 𝑇 𝑾𝑢 𝒖 𝑘 − 𝒖𝑓 𝑘=0 リカチ行列式を解き、次の定常ゲイン行列Kを決定する 𝑲 = 𝑾𝑢 + 𝑸𝑻 𝑯𝑸 −1 𝑸𝑻 𝑯𝑷 このKによる状態ベクトルフィードバックは、 𝒖 𝑘 − 𝒖𝑓 = −𝐊 𝒙(𝑘) − 𝒙𝑓 多入出力プラントが０型で一定入力(uf )に対して定常出力 (yf )が存在する場合 𝒙𝑓 = lim 𝑧 − 1 𝑿(𝑧) = 𝑰 − 𝑷 𝑧→1 𝒚𝑓 = 𝑪𝒙𝑓 −1 𝑸𝒖 𝑓 66

7.5 多変数制御系（非干渉化） 67 状態ベクトルフィードバック： 𝒖 𝑘 − 𝒖𝑓 = −𝐊 𝒙(𝑘) − 𝒙𝑓 ここで入力信号w(k)を導入 𝒖 𝑘 = 𝒖𝑓 + 𝐊𝒖𝑓 − 𝐊𝒙 𝑘 = 𝑰𝑚 + 𝐊 𝑰𝑛 − 𝑷 −1 𝑸 𝒖𝑓 − 𝐊𝒙 𝑘 𝒖 𝑘 = 𝐅𝒘 𝑘 − 𝐊𝒙 𝑘 𝒘 𝑘 : 設定入力信号 入力信号w(k)は、その定常値wfが出力の定常値yfと一致する入力信号で、次のようなる 𝒚𝑓 = 𝑮 1 𝒖𝑓 ただし 𝑮(𝑧) = 𝑪 𝑧𝑰𝑛 − 𝑷 −1 𝑸 これらより、所要のフィルタFは次のように決定される 𝐅 = 𝑰𝑚 + 𝐊 𝑰𝑛 − 𝑷 −1 𝑸 ∙ 𝑮(1)−1 (7-32) また、入力w(k)とし、出力y(k)とする伝達関数Gf(z)は次のようになる 𝑮𝒇 (𝑧) = 𝑮(𝑧) ∙ 𝑭 (7-33)

7.5 多変数制御系（非干渉化） 例）3連振動系 入力 出力 出力 68 𝑛 = 6, 𝑚 = 2 入力 2入力2出力系 y1 (t) 対象システム： u1 (t) y2 (t) b1 m1 k1 b2 m2 y3 (t) u2 (t) k2 b3 m3 k3 k1=1.;k2=2.;k3=1.;b1=0.01;b2=0;b3=0.01;m1=m2=m3=1.0 0 −𝑘1 0 𝐴= 𝑘1 0 0 1 −𝑏1 0 𝑏1 0 0 0 0 1.0 0 0 0 𝐵= 0 0 0 0 0 2.0 0 𝑘1 0 − 𝑘1 + 𝑘2 0 𝑘2 0 𝑏1 1 − 𝑏1 + 𝑏2 0 𝑏2 0 0 0 𝑘2 0 − 𝑘2 + 𝑘3 𝐶= 1.0 0 0 0 0 𝑏2 1 − 𝑏2 + 𝑏3 0 0 0 0 0 1.2 0 0 （B,Cの数値は書籍表4-1 (p58）

7.5 多変数制御系（非干渉化） 69 ｗ＝0.1 u1--->y1, y2 Gain= [[ 1.12363 1.63764 0.59306 0.41807 0.31329 0.11736] [ 0.56432 0.23531 -0.58738 0.64909 0.79793 1.03578]] ;w= 0.1 dc gain= [[2.5 2. ] [1.8 2.4]] F gain= [[ 1.96698 0.0525 ] [-0.58465 1.54959]] u2--->y1, y2 DGC_no5/ p145v2

7.5 多変数制御系（非干渉化） 70 ｗ＝0.01 u1--->y1, y2 Gain= [[ 2.80202 2.50703 0.81748 0.43021 0.30338 0.07786] [ 0.8595 0.20358 -1.07651 1.55657 2.00353 1.43021]] ;w= 0.01 dc gain= [[2.5 2. ] [1.8 2.4]] F gain= [[ 3.65033 0.22713] [-0.89226 2.64898]] u2--->y1, y2 DGC_no5/ p145v2

7.5 多変数制御系（非干渉化） 71 page145 “高橋安人,ディジタル制御 ,岩波書店,1987

7.5 多変数制御系（ふろく） 以下のことを興味本位から行ってみた １． ２変数プロセスでRの非対角項に値を入れて応答の変化をみた 𝑹= 𝑝1 0 0 𝑝2 𝑝1 𝑹= 𝛽 𝛼 𝑝2 Rが対角形であるのは応答に含まれる等比を入力だけから決めることを意味する。たとえばu1→y1, u1→y2の 応答成分がともに等比p1をもって減衰するものと考える。 ２． ２変数プロセス（表7-2）でObserverを作ってみた Observer Feedback gainは式(4-30)を用いた 𝐹 𝑇 = 𝑰 𝑹 𝑹2 ⋯ 𝑹𝑛−2 𝑹𝑛−1 𝑮−1 𝑛 (4-30) （時系列形の伝達関数に対応する） ３． ２変数プロセスで非干渉化とLQI制御を比べてみた 72

7.5 多変数制御系 その１(Rの非対角項） 73 R Matrixの対角項以外に も値を入れてみた 干渉 （その値の根拠はない） 干渉 若干、干渉の様子が変化した （干渉が少なくなった？） noise noise step外乱 (K=20でv1, k=20でv2) V=1を入力 DGC_no5/ p142v1

7.5 多変数制御系 その2 (Observer） 74 2変数プロセスのObserver 𝐹𝑇 = 𝑰 𝑹 𝑹2 ⋯ 𝑹𝑛−2 𝑹𝑛−1 𝑮−1 𝑛 Observer Feedback gainは式 (4-30)を用いた DGC_no5/ p142obs

75 Observe出力 推定値(青)とその状態値 (赤)を重ねて表示 ⚫ 初期値に誤差を 与えている XH=np.array( [[4.], [2.], [-5.], [5.], [-10.], [-15.], [10.], [20.], [-9.], [8.], [0.], [0.]]) 状態変数は12個ある が、今回6個表示する 状態変数の推定値が X03以降振動している DGC_no5/ p142obs

7.5 76 多変数制御系 その2 (ObserverのR対角項 ） 2変数プロセスのObserver 𝐹𝑇 = 𝑰 𝑹 𝑹2 ⋯ 𝑹𝑛−2 𝑹𝑛−1 𝑮−1 𝑛 R Matrixの対角項以外に も値を入れてみた （その値の根拠はない） DGC_no5/ p142obs

Observe出力 推定値(青)とその状態値 (赤)を重ねて表示 ⚫ 初期値に誤差を 与えている XH=np.array( [[4.], [2.], [-5.], [5.], [-10.], [-15.], [10.], [20.], [-9.], [8.], [0.], [0.]]) 状態変数は12個ある が、今回6個表示する 状態変数の推定値がX03 以降の振動が収まる 77

7.5 多変数制御系 その3 (非干渉化v.s.LQI制御） ｗ＝0.1 LQI制御でも行ってみた u1--->y1, y2 Gain= [[ 1.35907 0.14036 2.22431 2.24222 0.89284 0.53568 0.37947 0.09854] [-0.34479 1.0379 0.82048 0.20498 0.04001 1.81252 1.77227 1.30946]] ;w= 0.1 78 u2--->y1, y2 例プラントでの結果を見る限り、 非干渉化はLQI制御でも実現可能である。対 角項優位の方式とほぼ同様な結果となった

7.5 多変数制御系 その3 (非干渉化v.s.LQI制御） ｗ＝0.01 LQI制御でも行ってみた u1--->y1, y2 Gain= [[ 2.73457 0.21236 3.71937 2.85618 1.05706 0.53449 0.33772 0.06787] [-0.65441 1.96118 1.13197 0.19392 0.48752 3.19646 3.09367 1.63202]] ;w= 0.01 79 u2--->y1, y2 プログラムは6.5節で利用した ものとほぼ同じものとなる

つづく ここまで７章です。次回は8章以降に進みます。 • 80 （期間が開くかもしれません） やっとここまで進みました。ここらまで来ると書籍にはプログラムリストはなく、数式の 記述と簡単な解説、そして例題の処理結果の掲載になります。この処理結果には必要な条 件＆パラメータと計算結果の値が適切に書かれているため大変助かります。たぶん著者自 らがプログラムを作られたのでキーポイントを押さえているのだなと感心した次第です。 • ここで気が付いたのは、計算結果数値の微妙な違いです。たぶん計算を行った年代、すな わち計算機＆計算アルゴリズムの違いによるのかなと思います。 • 7章では、Rに非対角項を入れたり、Observerを作成したり、非干渉化でLQIとの比較を したりで、興味本位の項目を入れてみました。但し、プログラムなど含めこれら内容の信 頼性は、私が作ったものなので疑問です。ご注意ください。 • ９，１０章は少々内容がオムニバスのようになってきたので、私のスライドは、たぶん8 章で最後かなと思っています （Pythonのprogramは下記のGitHubに格納してあります） https://github.com/m4881shimojo/YTakahashi_Digital-control/tree/main

• https://github.com/m4881shimojo/YTakahashi_Digital-control/tree/main

ダイナミックプログラミング1/4 81 Riccati式の誘導とダイナミックプログラミング 対象は1次系とする： (6-1) 𝑥 𝑘 + 1 = 𝑝𝑥 𝑘 + 𝑞𝑢(𝑘) 𝑁−1 問題：Jを最小にする入力u(k)を求める。 𝐽 = ෍ 𝑊𝑥 𝑥 2 𝑘 + 1 + 𝑤𝑢2 (𝑘) (6-2) 𝐾=0 ダイナミックプログラミング：時点kにおける最適操作量u(k)は、k以前のことには関係 なくｋから先の最終ステップＮに至るまでの経過を最適にする(Bellman(1957)) 最終段から順次前へ最適化を進めていく 第N段のコスト（評価）： 𝑓(𝑁) = 𝑊𝑥 𝑥 2 𝑁 + 𝑤𝑢2 𝑁 − 1 これを最小にするu(N-1)を求める 𝜕𝑓Τ𝜕 𝑢 𝑁 − 1 = 0 𝑢 𝑁−1 =− 𝑞𝑊𝑥 𝑝 𝑥(𝑁 − 1) 𝑞2 𝑊𝑥 + 𝑤 これを(a)に代入、f(N)の最小値F(1)が次のようになる 最終段から 行う (a) (b)

ダイナミックプログラミング2/4 82 (c) 𝐹(1) = min 𝑓(𝑁) = 𝑅(𝑁 − 1) ∙ 𝑥 2 (𝑁 − 1) 𝑢(𝑁−1) 𝑅 𝑁−1 = 𝑝2 𝑊𝑥 𝑝2 𝑊𝑥 2 𝑞2 −− 2 𝑞 𝑊𝑥 + 𝑤 次に最後からの2ステップを考えると、 (d) 最終段-1 𝐹 2 = min 𝑝2 𝑥 2 𝑁 − 1 − 𝑤𝑢2 𝑁 − 2 + 𝐹 1 𝑤𝑢2 𝑘 𝑢 𝑁−2 = min 𝑢 𝑁−2 𝑊𝑥 + 𝑅 𝑁 − 1 𝑥 2 𝑁 − 1 − 𝑤𝑢2 𝑁 − 2 ) (e) ここで｛ ｝の項を次のようにおく 𝐻 𝑁 − 1 = 𝑊𝑥 + 𝑅 𝑁 − 1 (6-3) すると(e)は次のようになる 𝐹 2 = min H 𝑁 − 1 𝑥 2 𝑁 − 1 − 𝑤𝑢2 𝑁 − 2 ) (f) 𝑢 𝑁−2 (a), (c)にも(f)と同じ形が該当する。すなわち(a)のWxを 𝐻 𝑁 = 𝑊𝑥 𝐹(1) = min H 𝑁 𝑥 2 𝑁 − 𝑤𝑢2 𝑁 − 1 ) 𝑢 𝑁−1 (6-4) (g)

ダイナミックプログラミング3/4 83 (g)に対する最適制御式が(b)、そのときの最小コストが(c)と(d)であった。これに ならって、(g)と同形の(f)について計算すると次の結果となる 最適制御： 最小コスト： 𝑞𝐻(𝑁 − 1)𝑝 𝑢 𝑁−2 =− 2 𝑥(𝑁 − 2) 𝑞 𝐻(𝑁 − 1) + 𝑤 (h) 最終段-1 についての 最適制御 𝐹(2) = 𝑅(𝑁 − 1) ∙ 𝑥 2 (𝑁 − 1) 𝑅 𝑁−2 = 𝑝2 𝐻(𝑁 𝑝2 𝐻2 (𝑁 − 1)𝑞2 − 1) − 2 𝑞 𝐻(𝑁 − 1) + 𝑤 (i) (b)~(d)でWxをH(N)とおけば(h),(i)と完全に対応する 𝐻 𝑁 = 𝑊𝑥 + 𝑅(𝑁 − 2) 上式のR(N-2)に(i) を 代入すると 𝐻 𝑁 − 2 = 𝑊𝑥 + 𝑝2 𝐻(𝑁 𝑝2 𝐻2 (𝑁 − 1)𝑞2 − 1) − 2 𝑞 𝐻(𝑁 − 1) + 𝑤 (6-5)

ダイナミックプログラミング4/4 84 以下同様に、N-2,N-3…と続けることができる。ここで、N,N-1,…を0,1,2…にお きかえると次の差分式となる. 𝐻 𝑘+1 = 𝑝2 𝐻 𝑝2 𝐻2 𝑘 𝑞2 𝑘 − 2 + 𝑊𝑥 𝑞 𝐻 𝑘 +𝑤 一般形での 最適制御 (6-6) 𝐻 0 = 𝑊𝑥 𝑘 = 0,1,2, ⋯ 1次のリカチ式(Riccati Equation) lim 𝐻 𝑘 = 𝐻 𝑘→∞ (6-7) その定常解を使うと、最適制御が下記のように与えられる 𝑢 𝑘 = −𝐺 ∙ 𝑥(𝑘) 𝑞𝐻𝑝 𝐺= 2 𝑞 𝐻 𝑘 +𝑤 (6-8) (6-9)

85