ディジタル制御（1987年）その4 v2 pythonで記述してみる（6.線形2乗最適制御）

2024.3.15 ディジタル制御（1987年）その4 v2 pythonで記述してみる（6.線形2乗最適制御） 下 条 誠 電気通信大学名誉教授 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/

6. 6.1 線形2乗最適制御 行列リカチ(Riccati)式 6.2 慣性体のLQ制御 5章のPID制御法は直感的であった。6章 ではもっと厳密に2次系の評価関数数値 を最小化する。LQ制御(Linear quadratic 6.3 6.4 I動作を含むLQ制御 観測器を含むLQI制御 6.5 ０型プラントのLQI制御 optimal control)について行う 4

6.1 行列リカチ式（連続形） 最適レギュレータ問題 𝑑 𝒙 𝑡 = 𝑨𝒙 𝑡 + 𝑩𝑢(𝑡) 状態方程式： 𝑑𝑡 𝒚 𝑡 = 𝑪𝒙 𝑡 このシステムを制御するときの評価として次の2乗形を使う ∞ 𝐽 = න (𝒙𝑻 𝑸 𝒙 + 𝒖𝑻 𝑹𝒖)dt 0 𝑸 = 𝑸𝑇 > 0, 𝑹 = 𝑹𝑇 > 0 𝑸, 𝐑 ← 𝑤𝑒𝑖𝑔ℎ𝑡(𝑔𝑖𝑣𝑒) このJを最小にする入力u(t)を求める問題→ J を最小にするのは Pを求めu(t)を算出する 𝒖 = −𝑹−1 𝑩𝑇 𝑷 𝒙 である ただし、Pは下記リッカチ方程式(Riccati)を満たす解である 𝑷 = 𝑷𝑇 > 0 （連続形ではP、離散形ではHと表記） Riccati : 𝑨𝑇 𝑷 + 𝑷𝑨 − 𝑷𝑩𝑹−1 𝑩𝑇 𝑷 + 𝑸 = 0 𝑨, 𝑩, 𝑹, 𝑸 ← (𝑔𝑖𝑣𝑒𝑛) 5

6.1 行列リカチ式（離散形） 状態量Feedback system 6 𝒚 𝑘 𝒙 𝑘 + 1 = 𝑷𝒙 𝑘 + 𝒒𝒖 𝑘 𝒚 𝑘 = 𝒄𝒙(𝑘) 𝒖 𝑘 = −𝑮𝒙(𝑘) 𝑮 𝒙(𝑘) 𝑮 = 𝒘 + 𝒒𝑻 𝑯𝒒 −1 𝑻 𝒒 𝑯𝑷 u(k)を求める問題 G(feedback gain)を求める問題 H: リッカチ式の解を求めること 離散形ではH、連続形ではPと表記 Riccati行列： 𝑯 = 𝑷𝑻 𝑯𝑷 − 𝑷𝑻 𝑯𝒒 𝒘 + 𝒒𝑻 𝑯𝒒 𝑯 ∈ ℝ𝑛×𝑛 を求める −1 𝑻 𝒒 𝑯𝑷 + 𝑠𝑦𝑠𝑡𝑒𝑚 𝑔𝑖𝑣𝑒𝑛 → 𝑷, 𝒒 𝑷 ∈ ℝ𝑛×𝑛 , 𝒒 ∈ ℝ𝑛×𝑚 , 𝑾𝒙 ∈ ℝ𝑛×𝑛 , 𝒘 ∈ ℝ𝑚×𝑚 𝑾𝒙 𝑤𝑒𝑖𝑔ℎ𝑡 𝑔𝑖𝑣𝑒 → 𝒘, 𝑾𝒙 𝑯 ∈ ℝ𝑛×𝑛 , 𝑮 ∈ ℝ𝑚×𝑛

6.1 行列リカチ式（離散形：解法） H: リッカチ式の解を求める 7 反復計算による解法 行列リカチ(Riccati)式を反復計算すると安定解へ収束する 反復計算 𝑻 𝑯 𝟎 = 𝑾𝑥 = 𝒄 𝒄 𝑘 ←𝑘+1 注 初期値 𝑯 𝑘 が収束するまで繰り返す (6-13) 𝑯 𝑘 + 1 ← 𝑷𝑻 𝑯 𝑘 𝑷 − 𝑷𝑻 𝑯 𝑘 𝒒 𝑤 + 𝒒𝑻 𝑯 𝑘 𝒒 −1 𝑻 𝒒 𝑯 𝑘 𝑷 + 𝑾𝒙 𝑮 (6-14) 𝑯 → 𝑮 → 𝒖 𝑘 = −𝑮𝒙(𝑘) (6-15) ∞ 注) ここで評価関数Jの代表 的な形は離散系では 𝐽 = ෍ 𝑦 2 𝑘 + 1 + 𝑤𝑢2 𝑘 𝑦 2 𝑘 + 1 = 𝒙𝑇 𝒄𝑇 𝒄𝒙 → 𝑾𝒙 = 𝒄𝑇 𝒄 𝑘=0 DGC_no4/p101v1

6.1 行列リカチ式（解法例） 行列リカチ式の解法例を示す 例）解法例のシステム 𝑦(𝑡) 運動方程式： u(t) 𝑦ሷ 𝑡 + 𝑦(𝑡) ሶ =𝑢 𝑡 b ｚ変換式： m 𝐺𝑝 𝑧 = 𝐾𝑝 m=1, b=1 （5.3 状態方程式： 𝒙 𝑘 + 1 = 𝑷𝒙 𝑘 + 𝒒𝒖 𝑘 𝑷= 1 0 1−𝑝 𝑝 𝒒= 追従系の設計を参照） 𝒚 𝑘 = 𝒄𝒙(𝑘) 𝑝+𝑇 1−𝑝 𝑇:sampling interval , 𝒄= 1 𝑝 = 𝑒 −𝑇 ex) T=0.1, 𝑝 = 0.904837 ディジタル制御その３v2 5.3 追従系の設計を参照 𝑧+𝑞 𝑧−1 𝑧−𝑝 0 8

6.1 行列リカチ式(解の収束状況）𝑛 = 2, 𝑚 = 1 前頁システムのRiccati式を解きゲインGを求めた 𝑮 = 𝒘 + 𝒒𝑻 𝑯𝒒 w=0.1 9 −1 𝑻 𝒒 𝑯𝑷 ゲイン行列G(K1,K2) の収束状況 重み係数wを”0.1”と“1”の 二通りに変えた場合を示す 解Hの収束に従って Gも定まってくる w=1 DGC_no4/p103

6.1 行列リカチ式(解の収束状況） 10 書籍の図 “高橋安人,ディジタル制御 ,岩波書店,1987

6.1 行列リカチ式(解の収束状況）𝑛 = 2, 𝑚 = 1 11 同様に、重み係数wを”0.001”と“0.0001” に変えた場合を示す w=0.001 P matrix [[1. 0.09516258] [0. 0.90483742]] Q matrix [[0.00483742] （図6-3の値と同じ） [0.09516258]] 収束Num= 26 Gain= [[22.32400067 5.86734987]]] w=0.0001 P matrix [[1. 0.09516258] [0. 0.90483742]] Q matrix [[0.00483742] （図6-3の値と同じ） [0.09516258]] 収束Num= 14 Gain= [[52.46825213 9.47046297]] DGC_no4/p101v1

6.1 行列リカチ式 12 書籍のリスト page101 “高橋安人,ディジタル制御 ,岩波書店,1987

6.2 慣性体のLQ制御 6.1 行列リカチ式 6.2 慣性体のLQ制御 制御対象はtype 1型 6.3 制御対象はtype 0型 6.4 I動作を含むLQ制御 観測器を含むLQI制御 6.5 ０型プラントのLQI制御 1. 制御対象の運動方程式などを立て 2. 状態方程式（A, b)を記述し 3. 差分方程式の（P, q)に変換して 4. リカチ式を解き、ゲイン（Ｇ）を決め、 5. 対象システムに状態量をフィードバックする 13

6.2 慣性体のLQ制御 例 6.2節での手順をまとめてみた step1: 運動方程式： 14 𝑦(𝑡) u(t) b 𝑦ሷ 𝑡 + 𝑦(𝑡) ሶ =𝑢 𝑡 m システムの運動方程式 などの記述 step2: 状態方程式： ሶ 𝒙(𝑡) = 𝑨𝒙 𝑡 + 𝒃𝑢(𝑡) 差分式： 𝒙 𝑘 + 1 = 𝑷𝒙 𝑘 + 𝒒𝑢(𝑘) 𝑥1 𝑘 + 1 𝑥2 𝑘 + 1 = 1 1 − 𝑝 𝑥1 𝑘 0 𝑝 𝑥2 𝑘 状態方程式を求める step3: 差分方程式を求める （離散化） DGC_no4/p57_p101 step4: リカチ式： Riccati方程式を解く DGC_no4/p101v1 step5: 制御入力u(k)を決める step6: 制御出力を計算する 制御入力： 𝑯 = 𝑷𝑻 𝑯𝑷 − 𝑷𝑻 𝑯𝒒 𝒘 + 𝒒𝑻 𝑯𝒒 𝑠𝑦𝑠𝑡𝑒𝑚 𝑔𝑖𝑣𝑒𝑛 → 𝑷, 𝒒 𝒖 𝑘 = −𝑮𝒙(𝑘) −1 𝑻 𝒒 𝑯𝑷 + 𝑾𝒙 𝑤𝑒𝑖𝑔ℎ𝑡 𝑔𝑖𝑣𝑒 → 𝒘, 𝑾𝒙 𝑮 = 𝒘 + 𝒒𝑻 𝑯𝒒 −1 𝑻 𝒒 𝑯𝑷 DGC_no4/p101v1 応答： DGC_no4/p104 𝒚 𝑘 = 𝑪𝒙 𝑘 𝒙 𝑘 + 1 = 𝑷𝒙 𝑘 + 𝒒𝑢(𝑘) m=1, b=1 + 𝑝+𝑇−1 𝑢(𝑘) 1−𝑝

6.2 step1: 慣性体のLQ制御（対象とするシステム） 対象とするシステム 運動方程式： 伝達関数： （ｚ変換） 𝑦(𝑡) 15 u(t) b 𝑦ሷ 𝑡 + 𝑦(𝑡) ሶ =𝑢 𝑡 𝐺𝑝 𝑠 = 𝑌(𝑠) 1 = 𝑈(𝑠) 𝑠(𝑠 + 1) 𝐺𝑝 𝑧 = 𝐾𝑝 𝑝 = 𝑒 −𝑇 T：サンプリング間隔 𝑧+𝑞 𝑧−1 𝑧−𝑝 𝐾𝑝 = 𝑝 + 𝑇 − 1 𝑞= 1 − 𝑝 − 𝑝𝑇 𝑝+𝑇−1 m m: ｂ: y‘: u: m=1, b=1 質量 減衰係数 直線変位 力

6.2 step2-3: 慣性体のLQ制御（P,qを求める） 16 状態量 𝑷, 𝒒 を求める方法（２つ示す） 𝒖 𝑘 = −𝑮𝒙(𝑘) 𝒚 𝑘 𝒙 𝑘 + 1 = 𝑷𝒙 𝑘 + 𝒒𝒖 𝑘 𝒚 𝑘 = 𝒄𝒙(𝑘) 𝑮 𝒙(𝑘) 状態量をフィードバックする制御 方法１．連続形の状態方程式から求める 𝑑 𝒙 𝑡 = 𝑨𝒙 𝑡 + 𝒃𝑢(𝑡) 𝑑𝑡 方法２．伝達関数から求める 𝐺𝑝 𝑠 = 𝑌(𝑠) 1 = 𝑈(𝑠) 𝑠(𝑠 + 1) 𝑮 = 𝒘 + 𝒒𝑻 𝑯 𝑘 𝒒 状態量が必要 −1 𝑻 𝒒 𝑯 𝑘 𝑷 𝒙 𝑘 + 1 = 𝑷𝒙 𝑘 + 𝒒𝑢(𝑘) ディジタル制御(その２) 4.1節 𝒙 𝑘 + 1 = 𝑷𝒙 𝑘 + 𝒒𝑢(𝑘) 6.2節で述べる 𝒙 𝑘 + 1 = 𝑷𝒙 𝑘 + 𝒒𝑢(𝑘)

6.2 step2-3: 慣性体のLQ制御(方法1） 17 連続形の状態方程式から求める（ディジタル制御(その２) 4.1節） （１）1 入力系 A matrix [[ 0.00000 1.00000] [ 0.00000 -1.00000]] B matrix [[0] [1]] C matrix [[1.00000 0.00000]] P matrix [[1.00000000 0.09516258] [0.00000000 0.90483742]] 𝒙 𝑘 + 1 = 𝑷𝒙 𝑘 + 𝒒𝑢 𝑘 Q matrix [[0.00483742] [0.09516258]] 𝒚 𝑘 = 𝑐𝒙(𝑘) 𝑷 ∈ ℝ2×2 , 𝒒 ∈ ℝ2×1 , 𝒄 ∈ ℝ1×2 以下は付録です （2）2 入力系としてみた A matrix [[ 0.00000 1.00000] [ 0.00000 -1.00000]] B matrix [[1.00000 0.00000] [0.00000 1.00000]] P matrix [[1.00000000 0.09516258] [0.00000000 0.90483742]] Q matrix [[0.10000000 0.00483742] [0.00000000 0.09516258]] C matrix [[1.00000 0.00000] [0.00000 1.00000]] 𝑑 𝒙 𝑡 = 𝑨𝒙 𝑡 + 𝑩𝑢(𝑡) 𝑑𝑡 𝑷 ∈ ℝ2×2 , 𝒒 ∈ ℝ2×2 , 𝒄 ∈ ℝ2×2 𝑞11 𝒒𝒖 𝑘 = 𝑞 21 𝒚 𝑘 = 1 0 𝑞12 𝑢1 (𝑘) 𝑞22 𝑢2 (𝑘) 0 𝑦1 (𝑘) 1 𝑦2 (𝑘) となる？ 𝒙 𝑘 + 1 = 𝑷𝒙 𝑘 + 𝒒𝑢(𝑘) DGC_no4/p57_p101

6.2 慣性体のLQ制御(方法2） step2-3: 2. 伝達関数から求める（new） 𝑦(𝑡) u(t) b m 伝達関数を 変形して 𝑑 2 𝑦 𝑑𝑦 + =𝑢 𝑑𝑡 2 𝑑𝑡 𝐺𝑝 𝑠 = 𝑌(𝑠) 1 = 𝑈(𝑠) 𝑠(𝑠 + 1) ディジタル制御(その３) 5.3節 ｚ変換式： 𝑝 = 𝑒 −𝑇 𝑥1 = 𝑦 𝑥2 = 𝑦ሶ 𝑛 = 2, 𝑚 = 1 ラプラス変換式： 運動方程式： m=1,b=1 18 𝐺𝑝 𝑧 = 𝐾𝑝 𝑧+𝑞 𝑧−1 𝑧−𝑝 𝐾𝑝 = 𝑝 + 𝑇 − 1 𝑞= 1 − 𝑝 − 𝑝𝑇 𝑝+𝑇−1 状態量x1,x2を求める Gp(z)を直列結合で表す 𝐺𝑝 𝑧 𝐺𝑝 𝑧 1 𝑠 1+𝑠 𝑢(𝑘) hold 1 1+𝑠 𝐺𝑣 𝑧 𝑦(𝑘) 𝐺𝑣 𝑧 𝑥2 𝑘 1−𝑝 𝑧−𝑝 U(𝑧) 𝐺𝑣 𝑧 Τ𝐺𝑝 𝑧 𝑋2 𝑧 𝑝+𝑇−1 𝛼= 1−𝑝 𝛼𝑧 + 𝛽 𝑧−1 𝛽= Y(𝑧) 1 − 𝑝 − 𝑝𝑇 1−𝑝

6.2 慣性体のLQ制御(方法2） 19 𝑛 = 2, 𝑚 = 1 step2-3: 2. 伝達関数から求める（new） ①𝐺𝑣 𝑧 Τ𝐺𝑝 𝑧 に対して 𝛼𝑧 + 𝛽 𝑋1 𝑧 = 𝑋 𝑧 𝑧−1 2 ② 𝐺𝑣 𝑧 𝑧𝑋1 𝑧 = 𝑋1 𝑧 + 𝛼𝑧𝑋2 𝑧 + 𝛽𝑋2 𝑧 𝑌 𝑧 = 𝑋1 𝑧 に対して 𝑋2 𝑧 = 1−𝑝 𝑈 𝑧 𝑧−𝑝 𝑥1 𝑘 + 1 𝑥2 𝑘 + 1 上式より 𝑧𝑋2 𝑧 = 𝑝𝑋2 𝑧 + 1 − 𝑝 𝑈 𝑧 = 1 1 − 𝑝 𝑥1 𝑘 0 𝑝 𝑥2 𝑘 + 𝑝+𝑇−1 𝑢(𝑘) 1−𝑝 𝐺𝑝 𝑧 𝑷= 1 1−𝑝 0 𝑝 𝑝= 𝑒 −𝑇 𝒒= 𝑝+𝑇−1 1−𝑝 𝐺𝑣 𝑧 T：サンプリング間隔 U(𝑧) 1−𝑝 𝑧−𝑝 𝐺𝑣 𝑧 Τ𝐺𝑝 𝑧 𝑋2 𝑧 𝛼𝑧 + 𝛽 𝑧−1 Y(𝑧)

6.2 慣性体のLQ制御 step4-6: 20 𝑛 = 2, 𝑚 = 1 ステップ+ランプ応答 w=0.001 K1=22.324;K2=5.867 書籍と同じ値 𝑥1ሶ (𝑡) 0 = 𝑥2ሶ (𝑡) 0 1 𝑥1 𝑡 −1 𝑥2 𝑡 + 0 𝑢(𝑡) 1 DGC_no4/p104

6.2 慣性体のLQ制御 step4-6: 21 𝑛 = 2, 𝑚 = 1 ステップ+ランプ応答 w=0.0001 K1=52.467;K2=9.470 書籍と同じ値 𝑥1ሶ (𝑡) 0 = 𝑥2ሶ (𝑡) 0 1 𝑥1 𝑡 −1 𝑥2 𝑡 + 0 𝑢(𝑡) 1 DGC_no4/p104

6.2 慣性体のLQ制御 22 書籍の図 page104 目標入力が予め分 かっている場合、傾 斜入力より2刻み先 に操作量を上げた “高橋安人,ディジタル制御 ,岩波書店,1987

6.2 慣性体のLQ制御(2入力形） 23 以下の付録について検討 してみた これからの2力形については、書籍では言及がありません。私の興味本位での試みで 間違っている可能性大いにあります 以下は付録です （2）2 入力系としてみた A matrix [[ 0.00000 1.00000] [ 0.00000 -1.00000]] B matrix [[1.00000 0.00000] [0.00000 1.00000]] P matrix [[1.00000000 0.09516258] [0.00000000 0.90483742]] Q matrix [[0.10000000 0.00483742] [0.00000000 0.09516258]] C matrix [[1.00000 0.00000] [0.00000 1.00000]] 𝑑 𝒙 𝑡 = 𝑨𝒙 𝑡 + 𝑩𝑢(𝑡) 𝑑𝑡 𝑷 ∈ ℝ2×2 , 𝒒 ∈ ℝ2×2 , 𝒄 ∈ ℝ2×2 𝑞11 𝒒𝒖 𝑘 = 𝑞 21 𝒚 𝑘 = 1 0 𝑞12 𝑢1 (𝑘) 𝑞22 𝑢2 (𝑘) 0 𝑥1 (𝑘) 1 𝑥2 (𝑘) となる？ 𝒙 𝑘 + 1 = 𝑷𝒙 𝑘 + 𝒒𝑢(𝑘) DGC_no4/p57_p101

行列リカチ式(解の収束状況：2入力形） 24 ゲイン行列G(K11,K12,k21,k22)の収束状況 解Hの収束に従って Gも定まってくる 𝒙 𝑘 = 𝑥1 (𝑘) ∈ ℝ2×𝟏 𝑥2 (𝑘) u 𝑘 = 𝑢1 (𝑘) ∈ ℝ2×1 𝑢2 (𝑘) 𝑞11 𝒒= 𝑞 21 𝑮= 𝐾11 𝐾21 𝑞12 2×2 𝑞22 ∈ ℝ 𝐾12 ∈ ℝ2×2 𝐾22 DGC_no4/p101v3_m2

6.2 慣性体のLQ制御(応答の比較) （１）1 入力系として計算 𝒒 ∈ ℝ2×1 , 𝒘 ∈ ℝ1×1 DGC_no4/p104 25 （2）2 入力系として計算 重み w=0.01 𝒒 ∈ ℝ2×2 , 𝒘 ∈ ℝ2×2 𝒘= DGC_no4/p104_m2v2 𝑤 0 0 𝑤

6.2 慣性体のLQ制御(応答の比較) （１）1 入力系として計算 𝒒 ∈ ℝ2×1 , 𝒘 ∈ ℝ1×1 DGC_no4/p104 26 （2）2 入力系として計算 重み w=0.001 2 入力系としたほうが追従 は良さそうである 𝒒 ∈ ℝ2×2 , 𝒘 ∈ ℝ2×2 𝒘= DGC_no4/p104_m2v2 𝑤 0 0 𝑤

6.2 慣性体のLQ制御(m=1or2) 状態量Feedback 𝒚 𝑘 𝒙 𝑘 + 1 = 𝑷𝒙 𝑘 + 𝒒𝒖 𝑘 𝒚 𝑘 = 𝒄𝒙(𝑘) 𝒖 𝑘 = −𝑮𝒙(𝑘) 𝑮 𝒙(𝑘) 𝒒𝒖 𝑘 = −𝒒𝑮𝒙(𝑘) ∈ ℝ2×1 （m=1 or 2に関わらず、あたりまえか） （１）1 入力系 (m=1) 𝑮 = 𝐾1 , 𝐾2 ∈ ℝ1×2 𝒙 𝑘 = 𝑥1 (𝑘) ∈ ℝ2×𝟏 𝑥2 (𝑘) 𝒖 𝑘 = 𝑢(𝑘) ∈ ℝ1×1 𝑞1 𝒒 = 𝑞 ∈ ℝ2×1 2 （2）2 入力系 (m=2) 𝑥 (𝑘) 𝒙 𝑘 = 1 ∈ ℝ2×𝟏 𝑥2 (𝑘) u 𝑘 = 𝑢1 (𝑘) ∈ ℝ2×1 𝑢2 (𝑘) （この件で、1 or 2入力系という言い方が正しいかはわからない） 𝑮= 𝐾11 𝐾21 𝑞11 𝒒= 𝑞 21 𝐾12 ∈ ℝ2×2 𝐾22 𝑞12 2×2 𝑞22 ∈ ℝ 27

6.2 慣性体のLQ制御(m=1or2) 状態フィードバック：𝒖 𝑘 = −𝑮𝒙 𝑘 （１）1 入力系 (m=1) 𝒙 𝑘 = 𝑥1 (𝑘) ∈ ℝ2×𝟏 𝑥2 (𝑘) 𝒖 𝑘 = −𝑮𝒙 𝑘 = − 𝐾1 , 𝐾2 𝑞1 𝒒𝒖 𝑘 = − 𝑞 2 𝑮 = 𝐾1 , 𝐾2 ∈ ℝ1×𝟐 𝑥1 𝑘 𝑥2 𝑘 𝐾1 𝑥1 𝑘 +𝐾2 𝑥2 𝑘 = − 𝐾1 𝑥1 𝑘 +𝐾2 𝑥2 𝑘 𝑞1 𝒒 = 𝑞 ∈ ℝ2×1 2 ∈ ℝ1×1 𝑞1 𝐾1 𝑥1 𝑘 +𝐾2 𝑥2 𝑘 =− ∈ ℝ2×1 𝑞2 𝐾1 𝑥1 𝑘 +𝐾2 𝑥2 𝑘 （2）2 入力系 (m=2) 𝒙 𝑘 = 28 𝑥1 (𝑘) ∈ ℝ2×𝟏 𝑥2 (𝑘) 𝒖 𝑘 = −𝑮𝒙 𝑘 = − 𝐾11 𝐾21 𝑞11 𝑞12 𝒒𝒖 𝑘 = − 𝑞 21 𝑞22 𝑞 (𝐾 𝑥 = − 11 11 1 𝑞21 (𝐾11 𝑥1 𝑮= 𝐾11 𝐾21 𝐾12 𝑥1 𝑘 𝐾22 𝑥2 𝑘 𝐾12 ∈ ℝ2×𝟐 𝐾22 =− 𝑞11 𝒒= 𝑞 21 𝐾11 𝑥1 𝑘 +𝐾12 𝑥2 𝑘 𝐾21 𝑥1 𝑘 +𝐾22 𝑥2 𝑘 𝑞12 2×2 ∈ ℝ 𝑞22 ∈ ℝ2×1 𝐾11 𝑥1 𝑘 +𝐾12 𝑥2 𝑘 𝐾21 𝑥1 𝑘 +𝐾22 𝑥2 𝑘 𝑘 +𝐾12 𝑥2 𝑘 ) + 𝑞12 (𝐾21 𝑥1 𝑘 +𝐾22 𝑥2 𝑘 ) ∈ ℝ2×1 𝑘 +𝐾12 𝑥2 𝑘 ) + 𝑞22 (𝐾21 𝑥1 𝑘 +𝐾22 𝑥2 𝑘 )

6.3 I動作を含むLQ制御 6.1 行列リカチ式 6.2 慣性体のLQ制御 6.3 I動作を含むLQ制御 6.4 観測器を含むLQI制御 6.5 ０型プラントのLQI制御 制御対象はtype 1型 制御対象はtype 0型 29

6.3 I動作を含むLQ制御 30 I（積分）動作の導入 V(z) 外乱 応答 設定点 R(z) + E(z) ディジタル制御 Gc(z) U(z) + + サンプラホールド 付きプラントGp(z) Y(z) 図5-1 単一ループディジタル フィードバック制御系 1. 6.2節では、Type-1型の対象（z=1に極を持つ）の場合を行った 2. 6.3節では、Type-0型の対象について行う。 Type-0型では追跡 (Tracking)動作などで、オフセットを生 じる。これを除くため、制御に積分(I)動作を導入する

6.3 I動作を含むLQ制御 1. 6.2節：Type-1型の対象 例) 𝐺𝑝 𝑧 = 𝐾𝑝 𝑧+𝑞 𝑧−1 𝑧−𝑝 31 (z=1に極を持つ) 2. 6.３節：Type-0型の対象 例) Laplace変換式 𝐺𝑝 𝑠 = 1 今回のプラント 𝑇1 𝑠 + 1 𝑇2 𝑠 + 1 1 𝐺 𝑧 = 1 − 𝑧 −1 𝒵 ℒ −1 𝐺 𝑠 𝑠 ｚ変換式 𝑝1 = 𝑒 −𝑇Τ𝑇1 𝐺𝑝 𝑧 = 𝑝2 = 𝑒 −𝑇Τ𝑇2 𝑟 = 𝑇2 Τ𝑇1 ≠ 1 𝑏1 𝑧 + 𝑏2 𝑧 − 𝑝1 𝑧 − 𝑝2 𝑏1 = 1 − 𝑝1 − 𝑟𝑝2 Τ 1 − 𝑟 𝑏2 = 𝑝1 𝑝2 − 𝑝2 − 𝑟𝑝1 Τ 1 − 𝑟 次の頁からこの制御対象を例とする

6.3 I動作を含むLQ制御 Type-0型の対象 例) V(𝑧) 𝐺𝑝 𝑧 𝐺𝑝1 𝑧 𝑰: 𝑅(𝑧) + − 𝐸(𝑧) 𝑧 𝑧−1 𝐾0 + − 𝑋0 𝑧 + + 1 − 𝑝2 𝑧 − 𝑝2 U(𝑧) + + 32 𝐺𝑝2 𝑧 𝑋2 𝑧 𝑏1 𝑧 + 𝑏2 Τ 1 − 𝑝2 𝑧 − 𝑝1 𝐾2 𝐾1 コントローラ 図6-4 プロセスの状態ベクトルフィードバック制御 ① I動作に関する状態変数をx0とする ② Gp(z)を2つに分け、その出力を状態変数x1,x2とする 状態方程式： 𝑥0 𝑘 + 1 𝑥1 𝑘 + 1 𝑥2 𝑘 + 1 1 = 0 0 −𝑝1 𝑝1 0 −𝑝3 𝑝3 𝑝2 𝑥0 𝑘 𝑥1 𝑘 𝑥2 𝑘 −𝑏1 𝑟 𝑘+1 𝑏1 + 𝑢 𝑘 + 0 1 − 𝑝2 0 次頁に導出過程を示す Y(𝑧) 𝑋1 𝑧

6.3 I動作を含むLQ制御 𝑋0 𝑧 = 状態変数 の関係式 𝑧 𝑅 𝑧 − 𝑋1 𝑧 𝑧−1 𝑏1 𝑧 + 𝑏2 Τ 1 − 𝑝2 𝑋1 𝑧 = 𝑋2 𝑧 𝑧 − 𝑝1 𝑋2 𝑧 = 1 − 𝑝2 𝑈 𝑧 𝑧 − 𝑝2 𝑝1 = 𝑒 −𝑇Τ𝑇1 , 𝑝2 = 𝑒 −𝑇Τ𝑇2 𝑏1 = 1 − 𝑝1 − 𝑟𝑝2 Τ 1 − 𝑟 𝑏2 = 𝑝1 𝑝2 − 𝑝2 − 𝑟𝑝1 Τ 1 − 𝑟 𝑟 = 𝑇2 Τ𝑇1 ≠ 1 𝑧𝑋0 𝑧 = 𝑋0 𝑧 + 𝑝1 𝑋1 𝑧 + 𝑝3 𝑋2 𝑧 + 𝑏1 𝑈 𝑧 +𝑧𝑅 𝑧 展開整理 𝑧𝑋1 𝑧 =𝑝1 𝑋1 𝑧 + 𝑝3 𝑋2 𝑧 +𝑏1 𝑈 𝑧 𝑧𝑋2 𝑧 =𝑝2 𝑋2 𝑧 + 1 − 𝑝2 𝑈 𝑧 𝑝3 = 𝑏2 + 𝑏1 𝑝2 Τ 1 − 𝑝2 𝑥0 𝑘 + 1 = 𝑥0 𝑘 − 𝑝1 𝑥1 𝑘 − 𝑝3 𝑥2 𝑘 − 𝑏1 𝑢 𝑘 +𝑟 𝑘 + 1 差分式へ （離散化） 𝑥1 𝑘 + 1 = 𝑝1 𝑥1 𝑧 + 𝑝3 𝑥2 𝑘 +𝑏1 𝑢 𝑘 𝑥2 𝑘 + 1 = 𝑝2 𝑥2 𝑘 + 1 − 𝑝2 𝑢 𝑧 前頁の状態方程式へ 33

6.3 I動作を含むLQ制御（失敗） 34 （本書籍で述べた） これまでの方法では失敗する（→積分器Gain=0となる） 最適制御のためP, qを用いてRiccati式を適用、ゲインGを求める 𝑮 = 𝒘 + 𝒒𝑻 𝑯𝒒 −1 𝑻 𝒒 𝑯𝑷 #図 6-4 のシステム q1=-b1;q2=b1;q3=1-p2 状 態 方 程 式 q1=-0.09433, q2=0.09433, P=np.array([[1,-p1, -p3], q3=0.4866 [0,p1,p3], [0,0,p2]]) Q=np.array([[q1], [q2], [q3]]) C=np.array([[0,1,0]]) 積分器のゲインがゼロ 原 因 u(k)が有限値となるため、 評価関数Jが発散する 対 策 制御入力を、d(k)=u(k)-u(k-1) の1刻み変化とする DGC_no4/p108

6.3 I動作を含むLQ制御 (成功) u(k)の代わりにその一刻み変化 d(k)=u(k+1)-u(k) 35 を用いる 1. このためには、状態変化も一刻み変化を新たな状態変数とする 𝑥1 𝑘 + 1 = 𝑥1 𝑘 + 𝑠1 𝑘 + 1 = 𝑥1 𝑘 + 𝑝1 𝑠1 𝑘 + 𝑝3 𝑠2 𝑘 + 𝑏1 𝑑 𝑘 + 𝑒(𝑘) 𝑠1 𝑘 + 1 = 𝑥1 𝑘 + 1 − 𝑥1 𝑧 = 𝑥1 𝑘 + 𝑝1 𝑠1 𝑘 + 𝑝3 𝑠2 𝑘 + 𝑏1 𝑑 𝑘 + 𝑒 𝑘 = 𝑝1 𝑠1 𝑘 + 𝑝3 𝑠2 𝑘 + 𝑏1 𝑑 𝑘 + 𝑒(𝑘) 𝑠2 𝑘 + 1 = 𝑥2 𝑘 + 1 − 𝑥2 𝑧 = 𝑝2 𝑥2 𝑘 − 𝑥2 𝑘 − 1 = 𝑝2 𝑠2 𝑘 + 1 − 𝑝2 𝑑 𝑘 + 𝑒(𝑘) + 1 − 𝑝2 𝑢 𝑘 − 𝑢 𝑘 − 1 𝑒 𝑘 = 𝑣 𝑘 − 𝑣(𝑘 − 1) 外乱 2. また出力は以下のようにする 𝑦 𝑘 + 1 − 𝑦 𝑘 = 𝑥1 𝑘 + 1 − 𝑥1 𝑘 = 𝑠1 𝑘 + 1 𝑦 𝑘 = 𝑥1 𝑘 3. この結果 状態方程式： 𝑥1 𝑘 + 1 𝑠1 𝑘 + 1 𝑠2 𝑘 + 1 𝑦2 𝑘+1 = 1 = 0 0 𝑝1 𝑝1 0 𝒙𝑇 𝒄𝑇 𝒄𝒙 𝑝3 𝑝3 𝑝2 → 𝑥1 𝑘 𝑠1 𝑘 𝑠2 𝑘 𝑏1 𝑏1 + 1 − 𝑝2 1 𝑾𝒙 = 𝒄 𝒄 = 0 1 0 𝑇 𝑑 𝑘 + 𝑒(𝑘) 0 1 0 = 0 0 0 0 0 0 0 0

6.3 I動作を含むLQ制御 (成功) u(k)の代わりにその一刻み変化 を用いる d(k)=u(k+1)-u(k) 最適制御のためP,qを用いてRiccati式を適用、ゲインGを求める 36 𝑮 = 𝒘 + 𝒒𝑻 𝑯𝒒 −1 𝑻 𝒒 𝑯𝑷 #図 6-4 のシステム 積分器のゲイン q1=-0.09433, q1=-b1;q2=b1;q3=1-p2 状 P=np.array([[1,p1, p3], 態 [0,p1,p3], 方 [0,0,p2]]) 程 式 Q=np.array([[q2], q2=0.09433, q3=0.4866 [q2], [q3]]) C=np.array([[1,0,0]]) 積分器のゲインK0が求まった 𝒖 𝑘 = −𝑮𝒙 𝑘 = − 𝐾0 DGC_no4/p110v1 𝐾1 𝐾2 𝑥1 𝑘 𝑠1 𝑘 𝑠2 𝑘

6.3 I動作を含むLQ制御 (d(k)) 37 プロセスのLQI制御例(w=0.1) Direct feedback y(k) 書籍と同じ 1 [0.17810479] 2 [0.52294361] 3 [0.83124519] 4 [1.00603365] 5 [1.06048085] 6 [1.05002833] 7 [1.02256673] 8 [1.00222894] 9 [0.99397379] 10 [0.99385012] k=10以降は入 11 [1.06790114] 力r(t)を変え 12 [1.27961805] たので比較の 意味なし 13 [1.61338137] 14 [2.0160239] 15 [2.43997442]] DGC_no4/p111v2

6.3 I動作を含むLQ制御 (d(k)) 38 プロセスのLQI制御例(w=0.01) Direct feedback y(k) 1 [0.37521105] 2 [0.87584674] 3 [1.07066944] 4 [1.05124882] 5 [1.00266026] 6 [0.9888179] 7 [0.99483803] 8 [1.00052151] 9 [1.0014052] 10 [1.00043914] 11 [1.14993273] 12 [1.50026029] 13 [1.92866111] 14 [2.3492164] 15 [2.75027183] 書籍と同じ k=10以降は入 力r(t)を変え たので比較の 意味なし DGC_no4/p111v2

6.3 I動作を含むLQ制御（I動作） 応答結果（書籍） w=0.1 39 page111 w=0.01 Gain Gain 書籍と同じ結果となった “高橋安人,ディジタル制御 ,岩波書店,1987 出力数値のリストのみで、 グラフ表記は無し

6.4 観測器を含むLQI制御 6.1 行列リカチ式 6.2 慣性体のLQ制御 6.3 I動作を含むLQ制御 6.4 観測器を含むLQI制御 6.5 ０型プラントのLQI制御 （4.5 観測器による状態推定を参照のこと） 40

6.4 観測器を含むLQI制御 41 本節では「応答実験により求めた伝達関数」を用いたLQI制御例を示す 4.4 𝑅(𝑧) + 𝐸(𝑧) − 𝐺𝑝 𝑧 V(𝑧) 𝑰: 𝑧 𝑧−1 𝐾0 + − 𝑋0 𝑧 + 時系列からの伝達関数を参照 + 𝑏1 𝑧 + 𝑏2 𝑧 − 𝑝1 𝑧 − 𝑝2 U(𝑧) 𝑲 Y(𝑧) 𝒙(𝑧) コントローラ 𝐾𝑒 −𝑠𝐿 制御対象： 𝐺𝑝 𝑠 = 1 + 𝑇1 𝑠 1 + 𝑇2 𝑠 𝑏1 𝑧 5 + 𝑏2 𝑧 4 + 𝑏3 𝑧 3 +𝑏4 𝑧 2 + 𝑏5 𝑧 + 𝑏6 𝐺𝑝 𝑧 = 𝑧 6 + 𝑎1 𝑧 5 1. まず伝達関数Gp(z)を実験による応答データ(gi)で表現する ←(4.4節 時系列からの伝達関数) 2. 次に、伝達関数Gp(z)から状態方程式の表現とする 3. そして、LQI制御を行う （LQIのobserver例は6.5節で示す） また、4.5 観測器による状態推定を参照のこと

6.4 観測器を含むLQI制御 42 まず伝達関数Gp(z)を実験による応答データ(gi)で表現する 例） 𝑏1 𝑧 5 + 𝑏2 𝑧 4 + 𝑏3 𝑧 3 +𝑏4 𝑧 2 + 𝑏5 𝑧 + 𝑏6 𝐺𝑝 = 𝑧 6 + 𝑎1 𝑧 5 𝑎1 = −𝑝 𝑔1 , 𝑔2 , 𝑔3 , 𝑔4 , 𝑔5 , 𝑔6 (𝑛 = 6) 𝑏1 = 𝑔1, 𝑏2 = 𝑔2 − 𝑝𝑔1 , 𝑏3 = 𝑔3 − 𝑝𝑔2 , 𝑏4 = 𝑔4 − 𝑝𝑔3 , 𝑏5 = 𝑔5 − 𝑝𝑔4, 𝑏6 = 𝑔6 − 𝑝𝑔5 応答データによる伝達関数Gp(z)から、状態方程式形に変換する方法は２つある 方法１：可制御標準形（同伴形）の利用（3.3節） 可制御正準形 (同伴形)とする 0 0 0 𝑃𝑐 = 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 1 0 −𝑎1 0 0 0 𝑞𝑐 = 0 0 1 方法２：時系列からの伝達関数で示された形（4.4節） 0 0 0 𝑃= 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 1 0 −𝑎1 𝑔1 𝑔2 𝑔3 𝑞= 𝑔 4 𝑔5 𝑔6 𝐶𝑐 = 𝑏6 𝑏5 𝑏4 𝑏3 𝑏2 𝑏1 (3-14) 応答データからのP行列は ランク落ちで、Observer 構成不可？ 𝐶= 1 0 0 0 0 0 (4-20)

6.4 観測器を含むLQI制御(方法1) １．プロセス特性Gp(z)に時系列形(実験データ gi)を用いる (4.4節 43 時系列からの伝達関数) 実測した応答データ： 𝑔1 , 𝑔2 , 𝑔3 , 𝑔4 , 𝑔5 , 𝑔6 (𝑛 = 6) 𝑏1 𝑧 5 + 𝑏2 𝑧 4 + 𝑏3 𝑧 3 +𝑏4 𝑧 2 + 𝑏5 𝑧 + 𝑏6 𝐺𝑝 = 𝑧 6 + 𝑎1 𝑧 5 𝑎1 = −𝑝 (4-19) 𝑏1 = 𝑔1 , 𝑏2 = 𝑔2 − 𝑝𝑔1 , 𝑏3 = 𝑔3 − 𝑝𝑔2 , 𝑏4 = 𝑔4 − 𝑝𝑔3 , 𝑏5 = 𝑔5 − 𝑝𝑔4 , 𝑏6 = 𝑔6 − 𝑝𝑔5 ２．状態方程式を求める 0 0 0 𝑃𝑐 = 0 0 0 可制御正準形 (同伴形)とする 𝑥1 𝑥2 𝑥3 𝑥4 𝑥5 𝑥6 𝑘+1 𝑘+1 𝑘+1 𝑘+1 𝑘+1 𝑘+1 0 0 0 = 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 1 −𝑎1 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 1 0 −𝑎1 𝑥1 𝑥2 𝑥3 𝑥4 𝑥5 𝑥6 𝑘 𝑘 𝑘 𝑘 𝑘 𝑘 0 0 0 + u(k) 0 0 1 0 0 0 𝑞𝑐 = 0 0 1 𝐶𝑐 = 𝑏6 𝑦 𝑘 = 𝑏6 𝑏5 𝑏4 𝑏3 𝑏2 𝑏1 (3-14) 𝑏5 𝑏4 𝑏3 𝑏2 𝑏1 𝑥1 𝑥2 𝑥3 𝑥4 𝑥5 𝑥6 𝑘 𝑘 𝑘 𝑘 𝑘 𝑘

6.4 観測器を含むLQI制御(方法1) u(k)の代わりにその一刻み変化 d(k)=u(k+1)-u(k) 44 を用いる (6.3節と同様ゲインK0が求まらないため) 𝑥0 𝑘 + 1 = 𝑥0 𝑘 + 𝑟 𝑘 + 1 − y(k + 1) 𝑥1 𝑘 + 1 = 𝑥2 𝑘 𝑠1 𝑘 + 1 = 𝑥1 𝑘 + 1 − 𝑥1 𝑘 = 𝑥2 𝑘 − 𝑥2 𝑘 − 1 = 𝑠2 𝑘 𝑥2 𝑘 + 1 = 𝑥3 𝑘 𝑠2 𝑘 + 1 = 𝑠3 𝑘 𝑥3 𝑘 + 1 = 𝑥4 𝑘 𝑠3 𝑘 + 1 = 𝑠4 𝑘 𝑠4 𝑘 + 1 = 𝑠5 𝑘 𝑥4 𝑘 + 1 = 𝑥5 𝑘 𝑠5 𝑘 + 1 = 𝑠6 𝑘 𝑥6 𝑘 + 1 = −𝑎1 𝑥6 𝑘 + u(k) 𝑠6 𝑘 + 1 = −𝑎1 𝑥6 𝑘 − 𝑥6 𝑘 − 1 + u k − u k − 1 = −𝑎1 𝑠6 𝑘 + d(k) 出力を新しい状態変数で与えるため、次のような計算をする 𝑦 𝑘+1 −𝑦 𝑘 =𝑐 𝒙 𝑘+1 −𝒙 𝑘 = 𝑏6 𝑠1 𝑘 + 1 + 𝑏5 𝑠2 𝑘 + 1 + 𝑏4 𝑠3 𝑘 + 1 + 𝑏3 𝑠4 𝑘 + 1 + 𝑏2 𝑠5 𝑘 + 1 + 𝑏1 𝑠6 𝑘 + 1 = 𝑏6 𝑠2 𝑘 + 𝑏5 𝑠3 𝑘 + 𝑏4 𝑠4 𝑘 + 𝑏3 𝑠5 𝑘 + 𝑏2 𝑠6 𝑘 + 𝑏1 −𝑎1 𝑠6 𝑘 + d(k) 𝑦 𝑘 + 1 − 𝑦 𝑘 = 𝑏6 𝑠2 𝑘 + 𝑏5 𝑠3 𝑘 + 𝑏4 𝑠4 𝑘 + 𝑏3 𝑠5 𝑘 + 𝑏2 −𝑎1 𝑏1 𝑠6 𝑘 + 𝑏1 d(k) 次頁にまとめ

6.4 観測器を含むLQI制御(方法1) 45 この結果、次の状態方程式が得られる 𝑦 𝑘+1 𝑠1 𝑘 + 1 𝑠2 𝑘 + 1 𝑠3 𝑘 + 1 𝑠4 𝑘 + 1 𝑠5 𝑘 + 1 𝑠6 𝑘 + 1 𝑦 𝑘 = 1 0 1 0 0 = 0 0 0 0 0 0 0 観測器設計での重み ベクトルの計算 0 𝑏6 0 1 0 0 0 0 0 0 0 0 0 0 𝑏5 0 1 0 0 0 0 𝑏4 0 0 1 0 0 0 𝑏3 0 0 0 1 0 0 P 𝑦 𝑘 𝑠1 𝑘 𝑠2 𝑘 0 0 𝑠3 𝑘 𝑠4 𝑘 𝑠5 𝑘 𝑠6 𝑘 𝒄𝑷 𝒄𝑷2 重みベクトル f： 𝒇 = 𝐹(𝑷) 𝒄𝑷3 ⋮ 𝒄𝑷𝑛 2 𝑏2 −𝑎1 𝑏1 0 0 0 0 1 −𝑎1 𝑇 𝑇 𝑦 𝑘+1 =𝒙 𝒄 −1 0 0 ⋮ 0 1 𝑦 𝑠1 𝑠2 𝑠3 𝑠4 𝑠5 𝑠6 𝑘 𝑘 𝑘 𝑘 𝑘 𝑘 𝑘 𝑏1 0 0 + 0 𝑑(𝑘) 0 0 1 1 0 𝒄𝒙 → 0 𝑇 𝑾𝒙 = 𝒄 𝒄 = 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 応答実験により求めた伝達関数ではobserver構成不可？ observerについては、 4.5 観測器による状態推定を参照のこと 0 0 0 0 0 0 0

6.4 観測器を含むLQI制御(方法1) 46 step状外乱 (step:25で) a1=-p0 V=0.5を入力 b1=g1;b2=g2-p0*g1;b3=g3-p0*g2;b4=g4- （書籍と同じ） p0*g3;b5=g5-p0*g4;b6=g6-p0*g5 r(k)入力 # P=np.array([[1,0,b6,b5,b4,b3,b2-a1*b1], [0,0,1,0,0,0,0], [0,0,0,1,0,0,0], [0,0,0,0,1,0,0], [0,0,0,0,0,1,0], y(k) [0,0,0,0,0,0,1], [0,0,0,0,0,0,-a1]]) #7x7 Q=np.array([[b1], [0], [0], noise [0], [0], [0], [1]]) #7x1 # C=np.array([[1,0,0,0,0,0,0]]) #1x7 u(k) d(k) gainの値はk0､k1は同じだが、 それ以外は違う。試しに、書籍 と同じ値としたら発振した DGC_no4/p115v4

6.4 観測器を含むLQI制御(方法1) 47 ここでは指令値r(k)変化 への追従性を検査した step40でramp入力を 与えた step状外乱 (step:25で) V=0.5を入力 （書籍と同じ） noise DGC_no4/p115v4

6.4 観測器を含むLQI制御(方法1) 48 書籍の図 page115 DGC_no4/p115v4 の結果と異なる “高橋安人,ディジタル制御 ,岩波書店,1987

6.4 観測器を含むLQI制御(方法2) １．プロセス特性Gp(z)に時系列形(実験データ gi)を用いる (4.4節 49 時系列からの伝達関数) 実測した応答データ： 𝑔1 , 𝑔2 , 𝑔3 , 𝑔4 , 𝑔5 , 𝑔6 (𝑛 = 6) 𝑏1 𝑧 5 + 𝑏2 𝑧 4 + 𝑏3 𝑧 3 +𝑏4 𝑧 2 + 𝑏5 𝑧 + 𝑏6 𝐺𝑝 = 𝑧 6 + 𝑎1 𝑧 5 𝑎1 = −𝑝 (4-19) 𝑏1 = 𝑔1 , 𝑏2 = 𝑔2 − 𝑝𝑔1 , 𝑏3 = 𝑔3 − 𝑝𝑔2 , 𝑏4 = 𝑔4 − 𝑝𝑔3 , 𝑏5 = 𝑔5 − 𝑝𝑔4 , 𝑏6 = 𝑔6 − 𝑝𝑔5 ２．状態方程式を求める 0 0 0 𝑃= 0 0 0 𝑥1 𝑥2 𝑥3 𝑥4 𝑥5 𝑥6 𝑘+1 𝑘+1 𝑘+1 𝑘+1 𝑘+1 𝑘+1 0 0 0 = 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 −𝑎1 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 1 0 −𝑎1 𝑥1 𝑥2 𝑥3 𝑥4 𝑥5 𝑥6 𝑘 𝑘 𝑘 𝑘 𝑘 𝑘 𝑔1 𝑔2 𝑔3 𝑞= 𝑔 4 𝑔5 𝑔6 𝑔1 𝑔2 𝑔3 + 𝑔 u(k) 4 𝑔5 𝑔6 𝐶= 1 𝑦 𝑘 = 1 0 0 0 0 0 0 (4-20) 0 0 0 𝑥1 𝑥2 𝑥 0 3 𝑥4 𝑥5 𝑥6 𝑘 𝑘 𝑘 𝑘 𝑘 𝑘

6.4 観測器を含むLQI制御(方法2) 50 この結果次の状態方程式が得られる 𝑦 𝑘+1 𝑠1 𝑘 + 1 𝑠2 𝑘 + 1 𝑠3 𝑘 + 1 𝑠4 𝑘 + 1 𝑠5 𝑘 + 1 𝑠6 𝑘 + 1 1 0 0 = 0 0 0 0 𝑦 𝑘 = 1 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 そして、 𝑦 2 𝑘 + 1 = 𝒙𝑇 𝒄𝑇 𝒄𝒙 → 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 −𝑎1 𝑦 𝑘 𝑠1 𝑘 𝑠2 𝑘 𝑠3 𝑘 𝑠4 𝑘 𝑠5 𝑘 𝑠6 𝑘 𝑔1 𝑔1 𝑔2 + 𝑔3 𝑑(𝑘) 𝑔4 𝑔5 𝑔6 𝑦 𝑘 𝑠1 𝑘 𝑠2 𝑘 0 0 𝑠3 𝑘 𝑠4 𝑘 𝑠5 𝑘 𝑠6 𝑘 1 1 0 0 0 0 0 0 0 0 0 0 𝑾𝒙 = 𝒄𝑇 𝒄 = 0 1 0 0 0 0 0 0 = 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

6.4 観測器を含むLQI制御(方法2) 51 step状外乱 (step:25で) V=0.5を入力 （書籍と同じ） # P=np.array([[1,0,1,0,0,0,0], [0,0,1,0,0,0,0], [0,0,0,1,0,0,0], [0,0,0,0,1,0,0], [0,0,0,0,0,1,0], [0,0,0,0,0,0,1], [0,0,0,0,0,0,-a1]]) #7x7 Q=np.array([[0], [g1], noise [g2], [g3], [g4], [g5], [g6]]) #7x1 # C=np.array([[1,0,0,0,0,0,0]]) #1x7 DGC_no4/p115p72v1

6.5 ０型プラントのLQI制御 6.1 行列リカチ式 6.2 慣性体のLQ制御 6.3 I動作を含むLQ制御 6.4 観測器を含むLQI制御 6.5 ０型プラントのLQI制御 • 1入力１出力系(n=3,m=1) • 2入力2出力系(n=6,m=2) 52

6.5 ０型プラントのLQI制御 53 本節では、LQI制御を０型プラント一般に拡張する 𝑅(𝑧) + 𝐸(𝑧) − 𝐺𝑝 𝑧 V(𝑧) 𝑰: 𝑧 𝑧−1 𝐾0 + − 𝑋0 𝑧 + + 𝑏1 𝑧 + 𝑏2 𝑧 − 𝑝1 𝑧 − 𝑝2 U(𝑧) 𝑲 Y(𝑧) 𝒙(𝑧) コントローラ 1. 制御入力を一刻み差とする状態方程式の一般形を示す 2. 次の2例について例題を示す ① 1入力１出力系(n=3,m=1)の例 ② 2入力2出力系(n=6,m=2)の例 observerについては、 4.5 観測器による状態推定を参照のこと

6.5 ０型プラントのLQI制御(m=1) 1入出力系 一般形への拡張 評価関数： 54 ∞ ∑ で記すべきですが 𝐽 = න (𝑦 𝑻 𝑾 𝑦 + 𝑑 𝑻 𝑤d)dt 0 状態方程式： 𝒙 𝑘 + 1 = 𝑷𝒙 𝑘 + 𝒒𝑢(𝑘) 𝑦 𝑘 = 𝒄𝒙(𝑘) 𝒙 ∈ ℝ𝑛×1 𝑷 ∈ ℝ𝑛×𝑛 𝒒 ∈ ℝ𝑛×1 𝑢 ∈ ℝ 1刻み差： 𝑦∈ℝ 𝒄 ∈ ℝ1×𝑛 𝒔 𝑘+1 = 𝒙 𝑘+1 −𝒙 𝑘 =𝑷 𝒙 𝑘 −𝒙 𝑘−1 𝒔 𝑘 + 1 = 𝑷𝒔 𝑘 + 𝒒𝑑 𝑘 𝑦 𝑘 + 1 − 𝑦 𝑘 = 𝒄 𝒙 𝑘 + 1 − 𝒙(𝑘) 𝑦 𝑘 + 1 = 𝑦 𝑘 + 𝒄𝑷𝒔 𝑘 + 𝒄𝒒𝑑 𝑘 新たな 状態方程式： 𝑦 𝑘+1 𝒔 𝑘+1 = 1 0 𝒄𝑷 𝑦 𝑘 𝑷 𝒔 𝑘 𝑷1 ∈ ℝ 𝑛+1 × 𝑛+1 𝒄𝒒 + 𝒒 𝑑(𝑘) 𝒒1 ∈ ℝ 𝑛+1 ×1 +𝒒 𝑢 𝑘 −𝑢 𝑘−1 𝒔 ∈ ℝ𝑛×1

6.5 ０型プラントのLQI制御(m=1) 一般形への拡張 𝑦 𝑘+1 𝒔 𝑘+1 = 𝒄𝒒 + 𝒒 𝑑(𝑘) 𝒒1 𝒄𝑷 𝑦 𝑘 𝑷 𝒔 𝑘 1 0 𝑷1 行列リカチ式 𝑯 = 𝑷1 𝑻 𝑯𝑷1 − 𝑷1 𝑻 𝑯𝒒1 𝒘 + 𝒒1 𝑻 𝑯𝒒1 にP1, q1を適用： 𝑛+1 × 𝑛+1 𝑷1 ∈ ℝ 重み行列： 1 0 0 𝑾𝑥 = ⋮ 0 0 0 0 0 ⋮ 0 0 0 0 0 ⋮ 0 0 ⋯ ⋯ ⋯ ⋯ ⋯ ⋯ 𝑇 ゲイン： 𝑮 = 𝑤 + 𝒒1 𝑯𝒒1 𝑮 = 𝐾0 𝐾0 適用制御則： 𝑑 𝑘 = −𝐺 ⋯ 𝑦 𝑘 𝑠 𝑘 𝒒1 ∈ ℝ 0 0 0 ⋮ 0 0 0 0 0 ∈ ℝ 𝑛+1 ⋮ 0 0 −1 𝒒𝑇 𝑯𝑷 1 1 ∈ℝ −1 𝑛+1 ×1 55 𝑚=1 1入出力系 𝒒1 𝑻 𝑯𝑷1 + 𝑾𝒙 𝑤∈ℝ × 𝑛+1 𝑛+1 ×1 𝐾𝑛 = −𝐾0 y k − 𝐾1 𝑠1 − ⋯ − 𝐾𝑛 𝑠𝑛 𝑑∈ℝ

6.5 ０型プラントのLQI制御(例1) 対象システム： 𝑦 𝑘+1 𝒔 𝑘+1 振動系： 𝐺 𝑧 = 1 𝒄𝑷 𝑦 𝑘 = 0 𝑷 𝒔 𝑘 𝑷1 B 𝑧 𝑏0 𝑧 + 𝑏1 = 2 𝐴 𝑧 𝑧 + 𝑎1 𝑧 + 𝑎2 𝒄𝒒 + 𝒒 𝑑(𝑘) 𝒒1 -------------P, Q matrix-----------w= 0.1 P matrix [[ 1.0000000 0.0000000 0.0670000] 𝑷1 [ 0.0000000 0.0000000 1.0000000] [ 0.0000000 -0.9330000 1.8660000]] Q matrix [[0.0000000] 𝒒1 [0.0000000] [1.0000000]] C matrix [[1 0 0]] 書籍と同じ 収束Num= 36 Gain= [[ 1.8462060 -0.6149891 1.0270583]] #p118 振動系 a1=-1.866; a2=0.933;b0=0.0; b1=0.067 P=np.array([[0,1], [-a2,-a1]]) # Q=np.array([[0], [1]]) #2x1 C=np.array([[b1,b0]]) #1x2 #LQI 制御 P1=np.array([[1,-a1*b0,b1-a1*b0], [0,0,1], [0,-a2,-a1]]) #3x3 Q1=np.array([[b0], [0], [1]]) #3x1 F= [[ 1.0000000] [14.9253731] [27.8507463]] C1=np.array([[1,0,0]]) #1x3 観測器の重みベクトル f 56

6.5 ０型プラントのLQI制御(例1) 𝐺 𝑧 = Open-loop step応答 57 𝑛 = 3, 𝑚 = 1 B 𝑧 𝑏0 𝑧 + 𝑏1 = 2 𝐴 𝑧 𝑧 + 𝑎1 𝑧 + 𝑎2 左図には制御なしの step応答も示す 周期約25刻みの減衰振動である step状外乱 (step:25で) V=0.5を入力 （書籍と同じ） 操作量u(k)と制御量y(k)の位相差 noise が設定点近傍での振動を抑えるの に有効に役立っているのがわかる LQI-Gains [[ 1.8462060 -0.6149891 1.0270583]] 書籍と同じ DGC_no4/p118v1

6.5 ０型プラントのLQI制御(observer) 𝑛 = 3, 𝑚 = 1 58 観測器(observer)出力 （書籍になし） 推定値(青)とその状態値(赤)を重ねて表示 （4.5とは青赤逆） n=３のシステム 有限整定観測では３刻みで収束する f : feedback observerについては、4.5 観測器による状態推定を参照のこと DGC_no4/p118v1

6.5 ０型プラントのLQI制御(例1) 59 書籍の図 page118 [[ 1.8462060 -0.6149891 1.0270583]] “高橋安人,ディジタル制御 ,岩波書店,1987

6.5 ０型プラントのLQI制御(例2) 𝑛 = 6, 𝑚 =602 対象システム： y2 (t) y1 (t) u1 (t) b1 m1 k1 u2 (t) b2 m2 u3 (t) m3 k2 2入力2出力系 y3 (t) b3 k3 3連振動系の状態式 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 外力を質量m1とm3 に加えた 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 page 239 m1=m2=m3=1 2出力としてm1とm2 の変位にした （B,Cの数値は書籍のまま） 3連振動系でバネ定数を(k1=1,k2=2,k3=1)、減衰係数を(b1=0.01,b2=0,b3=0.01)とした

6.5 ０型プラントのLQI制御(例2) 61 𝑛 = 6, 𝑚 = 2 2入力2出力系 状態方程式： 𝒚 𝑘+1 𝒔 𝑘+1 = 𝑰 𝒄𝑷 𝒚 𝑘 𝟎 𝑷 𝒔 𝑘 + 𝒄𝑸 𝒅(𝑘) 𝑸 𝑷1 n：3連振動系の状態数(n=6) m：入出力数(m=2) 𝑸1 𝑷 ∈ ℝ𝑛×𝑛 , 𝑸 ∈ ℝ𝑛×𝑚 , 𝑪 ∈ ℝ𝑚×𝑛 𝑷𝟏 ∈ ℝ 𝑛+𝑚 × 𝑛+𝑚 , 𝑸𝟏 ∈ ℝ 𝑛+𝑚 ×𝑚 , 𝑮 ∈ ℝ𝑚×(𝑛+𝑚) 𝒚 ∈ ℝ𝑚×1 , 𝒔 ∈ ℝ𝑛×1 , 𝒅 ∈ ℝ𝑚×1 Riccati方程式： 𝑯 = 𝑷𝟏 𝑻 𝑯𝑷𝟏 − 𝑷𝟏 𝑻 𝑯𝑸𝟏 𝒘 + 𝑸𝟏 𝑻 𝑯𝑸𝟏 𝑯∈ℝ 𝑛+𝑚 × 𝑛+𝑚 , 𝑾𝒙 ∈ ℝ 1 0 0 0 𝑾𝒙 = 0 0 0 0 ex. ∞ 𝐽 = න (𝒚𝑻 𝑾𝒙 𝒚 + 𝒅𝑻 𝒘𝐝)dt 0 𝑛+𝑚 × 𝑛+𝑚 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 −1 𝑸𝟏 𝑻 𝑯𝑷𝟏 + 𝑾𝒙 , 𝒘 ∈ ℝ𝑚×𝑚 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ex. 𝒘= 𝑤11 0 0 𝑤22

6.5 ０型プラントのLQI制御(例2) 指令入力r(k)の与え方： X=np.zeros((n+m,1)) Xoff=np.array([[0],[0],[0],[0],[0],[0],[0],[0]]) 状態変数(n=6,m=2) 入力r(k)を状態変数Xのオフセットと考える for k in range(1,knum): # #d[k] 積分が入ってると、こうなるのかな？ #入力 R(k)をオフセットとして設定 Xoff[0,0]=rinp[k,0]; Xoff[1,0]=rinp[k,1] UD=-np.dot(G,(X-Xoff)) rinp[k,0]： k番目のr1指令値 rinp[k,1]： k番目のr2指令値 G:LQIのFB ゲイン d[k]=UD.T #計算のため置き換え if k==0:u[k]=d[k]# For Plot PLOTのための記録 else:u[k]=u[k-1]+d[k] #外乱入力はここに入れる？。d(k)に加算 if k==Noise1:VD[0,0]=1.0 #r1(k) これらの方法が 正しいかは不明 入力ノイズ （1刻み入力：stepを表現） else:VD[0,0]=0 if k==Noise2:VD[1,0]=1.0 #r2(k) else:VD[1,0]=0 X=np.dot(P1,X)+np.dot(Q1,(UD+VD)) 新たな状態量の算出 #XX0[k]=np.transpose(X) #for PLOT y[k]=np.transpose(np.dot(C1,X)) # PLOTのための記録 DGC_no4/p120 62

6.5 ０型プラントのLQI制御(例2) Tsample=0.5 #sampling time 63 𝑛 = 6, 𝑚 = 2 -------------P1 and Q1 matrix------------ k1=1.;k2=2.;k3=1.;b1=0.01;b2=0;b3=0.01 #三連振動系 #4.1 プラント微分方程式と等価の差分式を参照のこと A=np.array([[0,1.0,0.,0.,0.,0.], [-k1,-b1,k1,b1,0.,0.], [0.,0.,0.,1.,0.,0.], [k1,b1,-(k1+k2),-(b1+b2),k2,b2], [0.,0.,0.,0.,0.,1.], P1 matrix [[ 1. 0. [ 0. 1. [ 0. 0. [ 0. 0. [ 0. 0. [ 0. 0. [ 0. 0. [ 0. 0. 0.8805 0.4785 0.1142 0.0209 0.0053 0.0005] 0.1375 0.0251 0.7929 0.528 0.2637 0.0463] 0.8805 0.4785 0.1142 0.0209 0.0053 0.0005] -0.4576 0.8759 0.4169 0.1188 0.0402 0.0053] 0.1146 0.0209 0.6607 0.44 0.2198 0.0385] 0.4191 0.1188 -1.222 0.6565 0.7643 0.2194] 0.0049 0.0005 0.2198 0.0385 0.6581 0.4397] 0.038 0.0053 0.7643 0.2194 -1.242 0.6537]] [0.,0.,k2,b2,-(k2+k3),-(b2+b3)]]) B=np.array([[0.,0.], #p58 hyou4-1 [1.,0.], [0.,0.], [0.,0.], [0.,0.], [0.,2.]]) C=np.array([[1.,0.,0.,0.,0.,0.], [0.,0.,1.2,0.,0.,0.]]) P matrix [[ 0.8805 0.4785 0.1142 0.0209 0.0053 0.0005] [-0.4576 0.8759 0.4169 0.1188 0.0402 0.0053] [ 0.1146 0.0209 0.6607 0.44 0.2198 0.0385] [ 0.4191 0.1188 -1.222 0.6565 0.7643 0.2194] [ 0.0049 0.0005 0.2198 0.0385 0.6581 0.4397] [ 0.038 0.0053 0.7643 0.2194 -1.242 0.6537]] Q matrix [[0.1222 0.0001] [0.4785 0.0011] [0.0027 0.0099] [0.0209 0.0771] [0. 0.2345] [0.0005 0.8794]] Q1 matrix [[0.1222 0.0001] [0.0033 0.0119] [0.1222 0.0001] [0.4785 0.0011] [0.0027 0.0099] [0.0209 0.0771] [0. 0.2345] [0.0005 0.8794]] 収束Num= 55 G:LQI ゲイン Gain= [[ 1.3591 0.1404 2.2243 2.2422 0.8928 0.5357 0.3795 0.0985] [-0.3448 1.0379 0.8205 0.205 0.04 1.8125 1.7723 1.3095]] ;w= 0.1 DGC_no4/p120

6.5 ０型プラントのLQI制御(例2) 64 𝑛 = 6, 𝑚 = 2 2入力2出力系 1. 外力を質量m1（自由端）とm3に 加える2入力系とした 2. 2出力としてm1とm2の変位にした LQI-Gains G= [[ 1.3591 0.1404 2.2243 2.2422 0.8928 0.5357 0.3795 0.0985] noise noise [-0.3448 1.0379 0.8205 0.205 1.8125 1.7723 1.3095]] 0.04 このシステムでは一つの設定点が突 変した時、他の制御量への干渉が、 ほとんどないことがわかる k=45：r2(k)=1→2でy1(k)変化ほぼ無し DGC_no4/p120

6.5 ０型プラントのLQI制御(例2) 65 書籍の図 page120 “高橋安人,ディジタル制御 ,岩波書店,1987

6.5 ０型プラントのLQI制御(例2) 66 書籍の図 page119 LQI-Gains “高橋安人,ディジタル制御 ,岩波書店,1987

つづく ここまで6章です。次回は7章以降に進みます。 • 67 （期間が開くかもしれません） 今回Pythonで用いた変数などは、書籍とできるだけ同じ名前を用いました。しかし、章 が進み、以前に使用していたリスト組合せながら新たな内容を追記したところ、変数名 の重なりや、異なる宣言が組み合わされた結果、変な現象が発生してきました。 • このため、これまでのプログラムと内容の改良を行いました。今回Ver2はこの見直しの 結果です。 • また、書籍にプログラムリストのないものも多くなってきました。これらは記述リスト を組合せて私が作成しました。これらにも誤解やミスなどがありそうです。 （Pythonのprogramは下記のGitHubに格納してあります） https://github.com/m4881shimojo/YTakahashi_Digital-control/tree/main

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

68