ディジタル制御（1987年）その6 pythonで記述してみる（8.　同定および適応制御）

シリーズ終了：2024.5.12 ディジタル制御（1987年）その6 pythonで記述してみる（8. 同定および適応制御） シリーズ終了 下 条 誠 電気通信大学名誉教授 https://researchmap.jp/read0072509/ https://www.docswell.com/user/m_shimojo “高橋安人，ディジタル制御 ，岩波書店， 1987 https://github.com/m4881shimojo/YTakahashi_Digital-control/tree/main The University of Electro-Communications 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

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

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

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

8. 同定および適応制御 8.1 簡単な適応サーボ系 8.2 模型規範適応制御 8.3 MRASによる同定と適応制御への応用 8.4 模型追従適応制御 8.5 最小2乗同定アルゴリズム 4

8.1 簡単な適応サーボ系 往（無負荷） 1. Armが往復運動をする． Armの角 度(y)を０から１と変え，１で停止 して物体ｍを掴み，０の位置へ戻る 2. この時０→１は無負荷，１→０は物 体ｍの負荷が増える ０ １ ０ m:質量 3. 負荷によって，質量がｍから２ｍに なったとする 復（負荷共） この様にArmの力学が往と復により違う場合のサーボ系について検討する 5

8.1 簡単な適応サーボ系（モデル決定） 6 例：armが往復する環境で，往では物体無し，復では物体を把持する系を想定する 𝑚 𝑑 2 𝑦 ′ (𝑡) 𝑑𝑡 ′ 2 𝑑𝑦 ′ (𝑡) +𝑏 = 𝑢(𝑡) 𝑑𝑡 ′ 𝑡 = 𝑏Τ𝑚 𝑡 ′ 𝑦 = 𝑏 2 Τ𝑚 𝑦 ′ 𝑑2 𝑦 ′ 𝑑𝑦 ′ 往（無負荷）： 𝑚 ′ 2 + 𝑏 ′ = 𝑢 𝑑𝑡 𝑑𝑡 𝐺𝑝 𝑠 = 𝐺𝑝 𝑧 = 𝐾𝑝 1 𝑠 𝑠+1 m=1 b=1 m 𝑑2𝑦 ′ 𝑑𝑦 ′ 復（負荷共）： 2𝑚 ′ 2 + 𝑏 ′ = 𝑢 𝑑𝑡 𝑑𝑡 𝑎0 1 𝐺𝑝 𝑠 = 𝑎 = 0 𝑠 𝑠 + 𝑎0 2 𝑧+𝑞 𝑧−1 𝑧−𝑝 𝑝 = 𝑒 −𝑇 𝐾𝑝 = 𝑝 + 𝑇 − 1 b m: 質量 ｂ: 減衰係数 y‘: 直線変位 u: 力 u(t) 𝑦(𝑡) 𝐺𝑝 𝑧 = 𝐾𝑝 𝑧+𝑞 𝑧−1 𝑧−𝑝 𝑝 = 𝑒 −𝑎0 𝑇 𝑞= 1 − 𝑝 − 𝑝𝑇 𝑝+𝑇−1 𝐾𝑝 = −1 + 𝑝 + 𝑎0 𝑇 𝑎0 𝑞= 1 − 𝑝 − 𝑎0 𝑇𝑝 −1 + 𝑝 + 𝑎0 𝑇

Laplace式からｚ式への変換： 1 𝐺 𝑧 = 1 − 𝑧 −1 𝒵 ℒ −1 𝐺 𝑠 𝑠 (4-8) この式をGp(s)に対して実行する（１→２ → 3 → 4→ ５） 𝑎0 1 𝐺 𝑠 = 𝑝 Gp(s)： 𝑎 = 0 𝑠 𝑠 + 𝑎0 2 1. 部分分数に変換する： 1 𝑎0 1 1 𝑎0 1 𝐺𝑝 𝑠 = 2 = − + 2+ 𝑠 𝑠 𝑠 + 𝑎0 𝑎0 𝑠 𝑠 𝑠 + 𝑎0 2. 逆Laplace変換+z変換： 𝒵 ℒ −1 𝐺𝑝 𝑠 = 1 𝑠 1 𝑧 𝑎0 𝑇𝑧 𝑧 − + + 𝑎0 𝑧−1 𝑧−1 2 𝑧−𝑝 3. (1-z-1)変換： 𝐺𝑝 𝑧 = 1 − 𝑧 −1 1 𝑧 𝑎0 𝑇𝑧 𝑧 − + + 𝑎0 𝑧−1 𝑧−1 2 𝑧−𝑝 4. (1-z-1)変換の結果： 𝐺𝑝 𝑧 = 5. 結 果： 1 𝑎0 𝑇 𝑧−1 1 −1 + + = 𝑎0 𝑧−1 𝑧−𝑝 𝑎0 𝐺𝑝 𝑧 = 𝐾𝑝 𝐾𝑝 = , 𝑝 = 𝑒 −𝑎0𝑇 −1 + 𝑝 + 𝑎0 𝑇 𝑧 − 𝑝 − 𝑎0 𝑇𝑝 + 1 𝑧−1 𝑧−𝑝 𝑧+𝑞 𝑧−1 𝑧−𝑝 −1 + 𝑝 + 𝑎0 𝑇 𝑎0 𝑞= 1 − 𝑝 − 𝑎0 𝑇𝑝 −1 + 𝑝 + 𝑎0 𝑇 𝑝 = 𝑒 −𝑎0𝑇 7

8.1 簡単な適応サーボ系（パラメータ） （5.3 追従系の設計参照） Limiter R(z) 設定点 + E(z) 制御方式 Gc(z) 𝐺𝑐 𝑧 = 𝐾𝑐 U(z) 𝑧+𝑏 , 𝑧+𝑎 U’(z) out in 𝑎 <1 Robot Arm系 Gp(z) 𝐺𝑝 𝑧 = 𝐾𝑝 Y(z) 𝑧+𝑞 𝑧−1 𝑧−𝑝 1. パラメータは，Robot Arm系Gpと制御系Gcがある． ① Robot Arm系Gp: Kp, p, q ② 制御系Gc: Kc, a, b 2. これらは，それぞれ，無負荷と負荷で値が違う 3. Robot Arm系のパラメータは負荷によって決まる 4. よって，設計するパラメータは制御系のみになる．そこでLQサー ボの理論（ 6.2，7.2節）を参考に，操作量の動きが激しすぎず， のろまでもない適切な値を選ぶ 5. その結果，根（固有値）は，0.5±0.3j, 0 とした 8 応答

8.1 簡単な適応サーボ系（根の指定） 制御対象： 制御方式： 特性式： 𝐺𝑝 𝑧 = 𝐾𝑝 𝑧+𝑞 𝑧−1 𝑧−𝑝 (5-16) 𝐺𝑐 𝑧 = 𝐾𝑐 𝑧+𝑏 𝑧+𝑎 (5-17) 𝑎 <1 1 + 𝐺𝑝 𝑧 𝐺𝑐 𝑧 → 𝐹 𝑧 = 𝑧 3 + 𝑎 1 𝑧 2 + 𝑎 2 𝑧 + 𝑎(3) (5-18) 𝐹 𝑧 = (𝑧 − 𝑧1 )(𝑧 − 𝑧2 ) 𝑧 − 𝑧3 = 𝑧 3 − 𝑧1 + 𝑧2 + 𝑧3 + 𝑧1 𝑧2 + 𝑧2 𝑧3 + 𝑧3 𝑧1 𝑧 − 𝑧1 𝑧2 𝑧3 特性式の固有値(z1,z2,z3)を指定することで，制御式のa, b, Kc が以下の手順で決まる 𝑎 1 = a−p−1+k 𝑎 2 = −ap − a + p + K b + q a(1), a(2), a(3)は特性値の固有値 (z1, z2, z3) を指定すると決まる 𝑎 3 = ap + Kbq a(1),(2),a(3) 決定 𝐹 1 = 1 + 𝑎 1 + 𝑎 2 + 𝑎(3) 制御式のa, b, Kcが (K = 𝐾𝑝 𝐾c ) K 1 + b = 𝐹(1)/(1 + 𝑞) 1 a= 𝑎 1 𝑞 + 𝑎(3) + 𝑞 1 + 𝑝 − 𝐾 1 + 𝑏 𝑞 𝑝+𝑞 決まる K=a 1 +1+p−a 9

8.1 簡単な適応サーボ系 無負荷用のパラメータを一貫して使用 質量2倍 質量2倍 10 負荷時に負荷用の制御パラメータを使用 質量2倍 質量2倍 少し小細工をした ① 負荷時は，limiterの値を大きくした ② 負荷時は， Kcを書籍の値とした（計算値の2倍となる） DGC_no6/p150

8.1 簡単な適応サーボ系 11 ----------------無負荷Arm---------------p=0.90484,Kp=0.0048374,q=0.967 ----------------負荷Arm---------------p=0.95123,Kp=0.0024588,q=0.983 ----------------無負荷制御---------------Kc=108.87,b=-0.67182,a=0.378 ----------------負荷制御---------------Kc=451.29,b=-0.69104,a=0.396 “高橋安人,ディジタル制御 ,岩波書店,1987 page150

8.2 模型規範適応制御 8.1 簡単な適応サーボ系 8.2 模型規範適応制御 8.3 MRASによる同定と適応制御への応用 8.4 模型追従適応制御 8.5 最小2乗同定アルゴリズム MRASの概要説明と program List 12

8.2 模型規範適応制御 固 𝑢(𝑘) 定 𝑓(𝑘) 𝑦(𝑘) 系 線 − アルゴリズム ℎ(𝑘) − ො 𝑦(𝑘) 変 + 𝑒(𝑘) + 適応 可 𝑓 𝑘 − ℎ(𝑘) 適応ループ 非線形 𝑔(𝑘) (b) 超安定ループ 𝑔(𝑘) 系 (a) MRAS 形 13 受け身 回路 𝑔 𝑘 ℎ(𝑘) (c) 受け身性 h(𝑘) 模型規範適応系 (Model Reference Adaptive System) 𝑒 𝑘 = 𝑦 𝑘 − 𝑦(𝑘) ො 固定系 可変系 lim 𝑒 𝑘 = 0 𝑘→∞ 適応アルゴリズムにより 可変系を変えていく 8.2 同 定 プラント モデル 両出力の一致をもって，プラントパラメータとする 8.4 追 跡 モデル プラント プラントがモデル出力に近ずくように入力を加える model: software model

8.2 模型規範適応制御 このプロセスの同定を行う プロセスのモデル 14 𝐺 𝑧 = 𝑔1 𝑧 −1 + 𝑔2 𝑧 −2 + ⋯ + 𝑔𝑁 𝑧 −𝑁 (8-6) 𝐺෠ 𝑧, 𝑘 = 𝑔ො1 𝑘 𝑧 −1 + 𝑔ො2 𝑘 𝑧 −2 + ⋯ + 𝑔ො𝑁 𝑘 𝑧 −𝑁 (8-7) 𝑘 = 0,1,2, ⋯ 目 標 ①プロセスの応答 lim 𝑔ො𝑗 𝑘 = 𝑔𝑗 𝑘→∞ 𝑗 = 0,1,2, ⋯ , 𝑁 𝑦 𝑘 + 1 = 𝑔1 𝑢 𝑘 + 𝑔2 𝑢 𝑘 − 1 + ⋯ + 𝑔𝑁 𝑢 𝑘 − 𝑁 + 1 𝑁 (8-8) (8-9) = ෍ 𝑔𝑗 𝑢 𝑘 + 1 − 𝑗 𝐽=1 𝑁 ②モデルの応答 𝑦ො 𝑘 + 1 = ෍ 𝑔ො𝑗 (𝑘 + 1)𝑢 𝑘 + 1 − 𝑗 (8-10) 𝐽=1 𝑁 ③モデルの暫定応答 𝑦ො ° 𝑘 + 1 = ෍ 𝑔ො𝑗 (𝑘)𝑢 𝑘 + 1 − 𝑗 𝐽=1 ①，②，③を使って目標を達成する (8-11)

8.2 模型規範適応制御 2通りの出力差： 𝑒 𝑘 + 1 = 𝑦 𝑘 + 1 − 𝑦ො 𝑘 + 1 15 (8-12) 𝑒 ° (𝑘 + 1) = 𝑦 𝑘 + 1 - 𝑦ො ° 𝑘 + 1 その差： 𝑒 ° 𝑘 + 1 − 𝑒 𝑘 + 1 = 𝑦ො 𝑘 + 1 - 𝑦ො ° 𝑘 + 1 𝑁 𝑒 ° 𝑘 + 1 − 𝑒 𝑘 + 1 = ෍ 𝑔ො𝑗 𝑘 + 1 − 𝑔ො𝑗 𝑘 ∙𝑢 𝑘+1−𝑗 (8-13) 𝐽=1 適応ループの超安定性を保証する適応アルゴリズム 𝑔ො𝑗 𝑘 + 1 = 𝑔ො𝑗 𝑘 + 𝐾𝑗 𝑒(𝑘 + 1) ∙ 𝑢 𝑘 + 1 − 𝑗 𝑗 = 0,1,2, ⋯ , 𝑁 (8-14) 𝐾𝑗 > 0 ただし，e(k+1) は(8-12)より， 𝑔ො𝑗 𝑘 + 1 を用いているのでe(k+1)を 実測することができない． そこで，式(8-10)を用いてe(k+1)を算定する

8.2 模型規範適応制御 16 そこでモデルの暫定応答を用いてe(k+1)を求める (8-14)式まま 𝑔ො𝑗 𝑘 + 1 = 𝑔ො𝑗 𝑘 + 𝐾𝑗 𝑒(𝑘 + 1) ∙ 𝑢 𝑘 + 1 − 𝑗 𝑔ො𝑗 適応アルゴリズム 𝑁 (8-13)式に代入 𝑒 ° 𝑘 + 1 − 𝑒 𝑘 + 1 = ෍ 𝑔ො𝑗 𝑘 + 1 − 𝑔ො𝑗 𝑘 ∙𝑢 𝑘+1−𝑗 𝐽=1 𝑁 𝑒 ° 𝑘 + 1 − 𝑒 𝑘 + 1 = ෍ 𝐾𝑗 𝑒(𝑘 + 1) ∙ 𝑢 𝑘 + 1 − 𝑗 ∙𝑢 𝑘+1−𝑗 𝐽=1 𝑁 𝑒 ° 𝑘 + 1 − 𝑒 𝑘 + 1 = 𝑒 𝑘 + 1 ෍ 𝐾𝑗 𝑢2 𝑘 + 1 − 𝑗 𝐽=1 e(k+1)の算定 𝑒 𝑘+1 = 𝑒° 𝑘 + 1 2 1 + σ𝑁 𝐽=1 𝐾𝑗 𝑢 𝑘 + 1 − 𝑗 (8-15) 1. 推定値 𝑔ො𝑗 を式(8-14)のルールで調整 2. その推定値 𝑔ො𝑗 を用いてモデル出力推定値 𝑦ො が式(8-10)より決まる 3. この出力推定値 𝑦ො とプラント出力yとの差としてeが生じ，１．に戻る

8.2 模型規範適応制御 MRAS 同定テスト プログラムList page155 “高橋安人,ディジタル制御 ,岩波書店,1987 17

8.3 MRASによる同定と適応制御への応用 8.1 簡単な適応サーボ系 8.2 模型規範適応制御 8.3 MRASによる同定と適応制御への応用 8.4 模型追従適応制御 8.5 最小2乗同定アルゴリズム 前節のアルゴリズムを 用いて例題を行う 18

8.3 MRASによる同定と適応制御への応用(例１) 例１： MRASを使って，初期モデルを対象プラントにあわせる（同定） 対象プラント： プロセス特性Gp(z)に時系列形(実験データ gi)を用いる (𝑛 = 4) 𝑝 = 0.6, 𝑔1 = 0.1, 𝑔2 = 0.3, 𝑔3 = 0.48, 𝑔4 = 0.4 𝑝0 = 0, 𝑔01 = 0.5, 𝑔02 = 0.4, 𝑔03 = 0.3, 𝑔04 = 0.15 初期モデル： 状態方程式： 𝒙 𝑘 + 1 = 𝑷𝒙 𝑘 + 𝒒𝑢(𝑘) 0 0 𝑃= 0 0 1 0 0 0 0 1 0 0 0 0 1 𝑝 𝑔1 𝑔2 𝑞= 𝑔 3 𝑔4 𝐶= 1 0 0 0 固定系 可変系 同 定 プラント モデル 両出力の一致をもって，プラントパラメータとする 追 跡 モデル プラント プラントがモデル出力に近ずくように入力を加える 19

8.3 MRASによる同定と適応制御への応用(例１) 20 プラントgと同定前のモデルg0 g =np.array([0.1,0.3,0.48,0.4]) g0=np.array([0.5,0.4,0.3,0.15]) dcGain=1とするようにNormalizeした g = [0.0532 0.1596 0.2553 0.2128] g0= [0.2660 0.2128 0.1596 0.0798]] （書籍ではdcGain=10となっているmiss？） MRAS同定が終了後 g0= [0.0532 0.1596 0.2553 0.2128] p0_hat= 0.6000000742595987 MRAS Gainは同じ値とした 𝑘1 = 𝑘2 = 𝑘3 = 𝑘4 = 100 random入力を使っているので波形は毎回違う DGC_no6/p157v1

8.3 MRASによる同定と適応制御への応用(例１) 21 時系列ｇと時系列g0のステップ応答 MRAS同定が終了後 プラントの時系列ｇとモデルの 時系列g0は．ほぼ一致する （当たり前だが） DGC_no6/p157v1

8.3 MRASによる同定と適応制御への応用(例１) 22 page157 “高橋安人,ディジタル制御 ,岩波書店,1987

8.3 MRASによる同定と適応制御への応用(例2) 23 例2：I動作がある制御系で対象プラントを同定し状態フィードバック制御を行う 対象プラント： プロセス特性Gp(z)に時系列形(実験データ gi)を用いる 𝑝 = 0.6, 𝑔1 = 0.1, 𝑔2 = 0.3, 𝑔3 = 0.48, 𝑔4 = 0.4 𝑝0 = 0, 𝑔01 = 0.5, 𝑔02 = 0.4, 𝑔03 = 0.3, 𝑔04 = 0.15 初期モデル： (𝑛 = 4)

8.3 MRASによる同定と適応制御への応用(例2) 例2：I動作がある制御系で対象プラントを同定し状態フィードバック制御を行う step1：制御系を作る ① Plantの状態方程式を時系列形(gi)を用いて作る ② このPlant制御系にI動作を入れる ③ 状態フィードバック系を構成する step2：Observerを用いて状態フィードバックを行う ① PlantをMRAS を用いて同定する ② 同定モデルの状態量を用いて状態フィードバックを行う 24

8.3 MRASによる同定と適応制御への応用(例2) 25 step1：① Plantの状態方程式を時系列形(gi)を用いて作る 対象プラント： プロセス特性に時系列形(gi)を用いる 𝑝 = 0.6, 𝑥1 𝑘 + 1 𝑥2 𝑘 + 1 𝑥3 𝑘 + 1 𝑥4 𝑘 + 1 𝑔1 = 0.1, 𝑔2 = 0.3, 𝑔3 = 0.48, 𝑔4 = 0.4 𝑝 = −𝑎1 0 0 𝑃= 0 0 1 0 0 0 0 0 = 0 0 0 0 1 0 0 1 0 −𝑎1 1 0 0 0 0 1 0 0 𝑔1 𝑔2 𝑞= 𝑔 3 𝑔4 0 0 1 𝑝 𝑥1 𝑘 𝑥2 𝑘 𝑥3 𝑘 𝑥4 𝑘 𝑔1 𝑔2 + 𝑔 u(k) 3 𝑔4 𝐶= 1 0 0 𝑦 𝑘 = 1 0 (𝑛 = 4) 0 𝑥1 𝑘 𝑥 𝑘 0 0 2 𝑥3 𝑘 𝑥4 𝑘

8.3 MRASによる同定と適応制御への応用(例2) 26 step1：② このPlant制御系にI動作を入れる 6.3 I動作を含むLQ制御を参照 𝑑 𝑘 =𝑢 𝑘 −𝑢 𝑘−1 𝑦 𝑘+1 𝑠1 𝑘 + 1 𝑠2 𝑘 + 1 𝑠3 𝑘 + 1 𝑠4 𝑘 + 1 1 0 = 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 1 𝑝 𝑦 𝑘 𝑠1 𝑘 𝑠2 𝑘 𝑠3 𝑘 𝑠4 𝑘 状態FB 𝑔1 𝑔1 + 𝑔2 𝑑(𝑘) 𝑔3 𝑔4 𝑦 𝑘 = 1 0 𝑦 𝑘 𝑠1 𝑘 0 0 0 𝑠2 𝑘 𝑠3 𝑘 𝑠4 𝑘

式の導出 27 １．状態方程式 𝑥1 𝑘 + 1 = 𝑥2 𝑘 +𝑔1 𝑢(𝑘) 𝑥1 𝑘 + 1 𝑥2 𝑘 + 1 𝑥3 𝑘 + 1 𝑥4 𝑘 + 1 𝑥2 𝑘 + 1 = 𝑥3 𝑘 +𝑔2 𝑢(𝑘) 𝑥3 𝑘 + 1 = 𝑥4 𝑘 +𝑔3 𝑢(𝑘) 𝑥4 𝑘 + 1 = −𝑎1 𝑥4 𝑘 + 𝑔4 𝑢(𝑘) 0 0 = 0 0 1 0 0 0 0 0 1 0 0 1 0 −𝑎1 𝑥1 𝑘 𝑥2 𝑘 𝑥3 𝑘 𝑥4 𝑘 𝑔1 𝑔2 + 𝑔 u(k) 3 𝑔4 ２．差分に変換する 𝑠1 𝑘 + 1 = 𝑥1 𝑘 + 1 − 𝑥1 𝑘 = 𝑥2 𝑘 + 𝑔1 𝑢 𝑘 − 𝑥2 𝑘 − 1 − 𝑔1 𝑢 𝑘 − 1 = 𝑥2 𝑘 − 𝑥2 𝑘 − 1 + 𝑔1 𝑢 𝑘 − 𝑢 𝑘 − 1 = 𝑠2 𝑘 + 𝑔1 𝑑 𝑘 𝑠2 𝑘 + 1 = 𝑥2 𝑘 + 1 − 𝑥2 𝑘 = 𝑥3 𝑘 + 𝑔2 𝑢 𝑘 − 𝑥3 𝑘 − 1 − 𝑔2 𝑢 𝑘 − 1 = 𝑠3 𝑘 + 𝑔2 𝑑 𝑘 𝑠3 𝑘 + 1 = 𝑥3 𝑘 + 1 − 𝑥3 𝑘 = 𝑥4 𝑘 + 𝑔3 𝑢 𝑘 − 𝑥4 𝑘 − 1 − 𝑔3 𝑢 𝑘 − 1 = 𝑠3 𝑘 + 𝑔3 𝑑 𝑘 𝑠4 𝑘 + 1 = 𝑥4 𝑘 + 1 − 𝑥4 𝑘 = −𝑎1 𝑥4 𝑘 + 𝑔4 𝑢 𝑘 −𝑎1 𝑥4 𝑘 − 1 − 𝑔4 𝑢 𝑘 − 1 = −𝑎1 𝑠4 𝑘 + 𝑔4 𝑑 𝑘 ３．出力を新しい状態変数で与える 𝑦 𝑘+1 −𝑦 𝑘 =𝑐 𝒙 𝑘+1 −𝒙 𝑘 = 𝑥1 𝑘 + 1 − 𝑥1 𝑘 = 𝑠1 𝑘 + 1 = 𝑠2 𝑘 + 𝑔1 𝑑 𝑘 この結果をまとめたのが前頁の状態方程式

8.3 MRASによる同定と適応制御への応用(例2) step1：③ 状態フィードバック系を構成する 状態FB 状態FB則より： 𝑑 𝑘 = −𝑮 𝑦 𝑘 𝒔 𝑘 = − 𝐾0 𝑦 𝑘 + 𝐾1 𝑠1 𝑘 + 𝐾2 𝑠2 𝑘 + 𝐾3 𝑠3 𝑘 + 𝐾4 𝑠4 𝑘 FBゲインGは，Riccati方程式より求める 𝑦 𝑘+1 𝑠1 𝑘 + 1 𝑠2 𝑘 + 1 𝑠3 𝑘 + 1 𝑠4 𝑘 + 1 1 0 = 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 1 𝑝 𝑦 𝑘 𝑠1 𝑘 𝑠2 𝑘 𝑠3 𝑘 𝑠4 𝑘 𝑔1 𝑔1 + 𝑔2 𝑑(𝑘) 𝑔3 𝑔4 𝑦 𝑘 = 1 0 0 𝑦 𝑘 𝑠1 𝑘 0 0 𝑠2 𝑘 𝑠3 𝑘 𝑠4 𝑘 28

8.3 MRASによる同定と適応制御への応用(例2) 29 １．差分入力d(k) 𝑦 𝑘 𝑑 𝑘 = −𝑮 = −𝐾0 𝑦 𝑘 − 𝐾1 𝑠1 𝑘 − 𝐾2 𝑠2 𝑘 − 𝐾3 𝑠3 𝑘 − 𝐾4 𝑠4 𝑘 𝒔 𝑘 状態FBの制御則 ２．この d(z) をZ領域で表すと 𝐷 𝑧 = 1 − 𝑧 −1 𝑈 𝑧 𝑆𝑖 𝑧 = 1 − 𝑧 −1 𝑋𝑖 𝑧 , 𝑖 = 1,2, ⋯ , 𝑛 ３．これを変形すると次のようになる 1 − 𝑧 −1 𝑈 𝑧 = 𝐷(𝑧) = −𝐾0 𝑌 𝑧 − 𝐾1 𝑆1 𝑧 − ⋯ − 𝐾𝑛 𝑆𝑛 𝑧 1 − 𝑧 −1 𝑈 𝑧 = −𝐾0 𝑌 𝑧 − 𝐾1 1 − 𝑧 −1 𝑋1 𝑧 ⋯ − 𝐾𝑛 1 − 𝑧 −1 𝑋𝑛 𝑧 𝑈 𝑧 = −𝐾0 (𝑛 = 4) ← 𝑆𝑖 𝑧 = 1 − 𝑧 −1 𝑋𝑖 𝑧 , ← 1 − 𝑧 −1 で割る 𝑧 𝑌 𝑧 − 𝐾1 𝑋1 𝑧 + ⋯ + 𝐾𝑛 𝑋𝑛 𝑧 𝑧−1 I(積分)動作 ４．設定点R(z)を付加すると次のようになる 𝑛 𝑧 𝑈 𝑧 = 𝐾0 𝑅 𝑧 − 𝑌(𝑧) − ෍ 𝐾𝑖 𝑋𝑖 𝑧 𝑧−1 𝑖=1 前頁の構成と一致する (6-26)

8.3 MRASによる同定と適応制御への応用(例2) 30 step2：① PlantをMRAS を用いて同定する 対象プラント： プロセス特性Gp(z)に時系列形(実験データ gi)を用いる 例1）と同じ 𝑝 = 0.6, 𝑔1 = 0.1, 𝑔2 = 0.3, 𝑔3 = 0.48, 𝑔4 = 0.4 𝑝0 = 0, 𝑔01 = 0.5, 𝑔02 = 0.4, 𝑔03 = 0.3, 𝑔04 = 0.15 初期モデル： (𝑛 = 4)

8.3 MRASによる同定と適応制御への応用(例2) 31 step2：② 同定モデルの状態量を用いて状態フィードバックを行う 入力のd(k)をu(k)に書き直す 定義より ： 𝑑 𝑘 =𝑢 𝑘 −𝑢 𝑘−1 状態FBより： 𝑑 𝑘 = − 𝐾0 𝑦 𝑘 + 𝐾1 𝑠1 𝑘 + 𝐾2 𝑠2 𝑘 + 𝐾3 𝑠3 𝑘 + 𝐾4 𝑠4 𝑘 𝑥0 𝑘 → 𝑦 𝑘 𝑠1 𝑘 = 𝑥1 𝑘 − 𝑥1 𝑘 − 1 𝑠2 𝑘 = 𝑥2 𝑘 − 𝑥2 𝑘 − 1 𝑠3 𝑘 = 𝑥3 𝑘 − 𝑥3 𝑘 − 1 𝑠4 𝑘 = 𝑥4 𝑘 − 𝑥4 𝑘 − 1 𝑦 𝑘+1 𝑠1 𝑘 + 1 𝑠2 𝑘 + 1 𝑠3 𝑘 + 1 𝑠4 𝑘 + 1 1 0 = 0 0 0 𝑛 𝑢 𝑘 = 𝐾0 𝑥0 𝑘 − ෍ 𝐾𝑖 𝑥𝑖 𝑘 𝑖=1 ただし，x0はI動作に含まれる状態変数で次のようになる 𝑥0 𝑘 = 𝑥0 𝑘 − 1 + 𝑟 𝑘 − 𝑦(𝑘) 𝑥0 𝑘 → 𝑦 𝑘 0 0 0 0 0 1 1 0 0 0 0 0 𝑦 𝑘 0 0 𝑠1 𝑘 0 0 𝑠2 𝑘 0 1 𝑠3 𝑘 0 𝑝 𝑠4 𝑘 𝑔1 𝑔1 + 𝑔2 𝑑(𝑘) 𝑔3 𝑔4

8.3 MRASによる同定と適応制御への応用(例2) step2：② 同定モデルの状態量を用いて状態フィードバックを行う 𝑰動作 𝑅(𝑧) + − 𝑧 𝑧−1 𝐾0 + U(𝑧) Y(𝑧) 𝐺𝑝 𝑧 − 𝑋0 𝑧 Observer 𝑛 − ෍ 𝐾𝑖 𝑥ො𝑖 𝑘 ෝ (𝑧) 𝒙 𝑲 𝑖=1 𝐾1 , 𝐾2 , 𝐾3 , 𝐾4 𝑛 Observer: MRASでの同定model 𝑧 𝑈 𝑧 = 𝐾0 𝑅 𝑧 − 𝑌(𝑧) − ෍ 𝐾𝑖 𝑋𝑖 𝑧 𝑧−1 𝑖=1 𝑢 𝑘 = 𝐾0 𝑥0 𝑘 − 𝐾1 𝑥ො1 𝑘 + 𝐾2 𝑥ො2 𝑘 + 𝐾3 𝑥ො3 𝑘 + 𝐾4 𝑥ො4 𝑘 𝑥0 𝑘 = 𝑥0 𝑘 − 1 + 𝑟 𝑘 − 𝑦(𝑘) 𝐺 = 𝐾0 , 𝐾1 , 𝐾2 , 𝐾3 , 𝐾4 𝑑 𝑘 = −𝑮 “G”は，Riccati方程式解より求める 𝑦 𝑘 𝒔 𝑘 32

8.3 MRASによる同定と適応制御への応用(例2) 33 適応動作は制御Gain 変化が一定範囲内の とき行うようにした (25～ ) TEST （TEST期間中は除く） DISTURBANCE 外乱は期間中Step状 の値を加えた (50～100) DGC_no6/p159

8.3 MRASによる同定と適応制御への応用(例2) 34 モデルパラメータの推移 対象プラント： 𝑝 = 0.6, 𝑔1 = 0.1, 𝑔2 = 0.3, 𝑔3 = 0.48, 𝑔4 = 0.4 初期モデル： 𝑝0 = 0, 𝑔01 = 0.5, 𝑔02 = 0.4, 𝑔03 = 0.3, 𝑔04 = 0.15 DGC_no6/p159

8.3 MRASによる同定と適応制御への応用(例2) 35 page159 “高橋安人,ディジタル制御 ,岩波書店,1987

8.4 模型追従適応制御 8.1 簡単な適応サーボ系 8.2 模型規範適応制御 8.3 MRASによる同定と適応制御への応用 8.4 模型追従適応制御 8.5 最小2乗同定アルゴリズム MRASには「同定」と 「追跡」がある この節は固定モデルの 挙動を追跡する話し 36

8.4 模型追従適応制御 37 摩擦係数が時間と共に変わるサーボ機構を， 有限整定応答の固定モデル挙動に追跡させる 制御例） b 運動方程式： 𝑦ሷ 𝑡 + 𝑦(𝑡) ሶ =𝑢 𝑡 伝達関数： 𝐺𝑝 𝑧 = 𝐾𝑝 𝑝 = 𝑒 −𝑇 u(t) 𝑦(𝑡) m m=1 b=1 𝑧+𝑞 𝑧−1 𝑧−𝑝 1 − 𝑝 − 𝑝𝑇 𝑞 = 𝐾𝑝 = 𝑝 + 𝑇 − 1 𝑝+𝑇−1 T：サンプリング間隔 m: ｂ: y‘: u: 質量 減衰係数 直線変位 力 固定系 可変系 同 定 プラント モデル 両出力の一致をもって，プラントパラメータとする 追 モデル プラント プラントがモデル出力に近ずくように入力を加える 跡

8.4 模型追従適応制御 38 𝑝 = 𝑒 −𝑏𝑇 𝑐0 = 𝑝 + 𝑏𝑇 − 1 Τ𝑏 2 𝑐1 = 1 − 𝑝 − 𝑝𝑏𝑇 Τ𝑏 2 𝑞0 = 𝑐0 Τ 1 − 𝑝 𝑞1 = 𝑐1 Τ 1 − 𝑝 特性式： 𝑧 2 + 𝐾1 1 − 𝑝 − 𝑝 − 1 + 𝐾0 𝑐0 𝑧 + 𝑝 − 𝐾1 1 − 𝑝 + 𝐾0 𝑐0 = 0 特性式有限整定→z2=0 より 閉ループ伝達関数： 𝐾0 = 𝑏 𝑇(1 − 𝑝) 1 − 𝑝 − 𝑏𝑝2 𝑇 𝐾1 = 𝑇(1 − 𝑝)2 （FBゲイン） 𝑌(𝑧) = 𝑧 −1 𝑏0 + 𝑏1 𝑧 −1 𝑅(𝑧) (8-22) （有限整定制御での） 𝑏0 = 𝑝 + 𝑏𝑇 − 1 𝑏𝑇(1 − 𝑝) 𝑏1 = 1 − 𝑝 − 𝑝𝑏𝑇 𝑏𝑇(1 − 𝑝)

前頁式の導出（Plant: P,Qを算出） 𝑞0 𝑧 + 𝑞1 𝑋 𝑧 = 𝑋 𝑧 ① X1に対して 1 𝑧−1 2 ② X2に対して 𝑋2 𝑧 = 1−𝑝 𝑈 𝑧 𝑧−𝑝 𝑧𝑋1 𝑧 = 𝑋1 𝑧 + 𝑞0 𝑧𝑋2 𝑧 + 𝑞1 𝑋2 𝑧 𝑥1 𝑘 + 1 = 𝑥1 𝑘 + 𝑞0 𝑥2 𝑘 + 1 + 𝑞1 𝑥2 𝑘 𝑧𝑋2 𝑧 = 𝑝𝑋2 𝑧 + 1 − 𝑝 𝑈 𝑧 𝑥2 𝑘 + 1 = 𝑝𝑥2 𝑘 + 1 − 𝑝 𝑢 𝑘 ①x1(k+1)に②x2(k+1)を代入して整理する 𝑥1 𝑘 + 1 = 𝑥1 𝑘 + 𝑞0 𝑝𝑥2 𝑘 + 1 − 𝑝 𝑢 𝑘 + 𝑞1 𝑥2 𝑘 = 𝑥1 𝑘 + 𝑞0 𝑝 + 𝑞1 𝑥2 𝑘 + 𝑞0 1 − 𝑝 𝑢 𝑘 𝑝𝑐0 𝑐1 𝑝 𝑝 + 𝑏𝑇 − 1 + 1 − 𝑝 − 𝑝𝑏𝑇 𝑝2 − 2𝑝 + 1 1 − 𝑝 𝑞0 𝑝 + 𝑞1 = + = = = 2 1−𝑝 1−𝑝 1 − 𝑝 𝑏2 1 − 𝑝 𝑏2 𝑏 𝑝 = 𝑒 −𝑏𝑇 𝑐0 𝑝 + 𝑏𝑇 − 1 𝑞0 1 − 𝑝 = 1−𝑝 = 1−𝑝 𝑏2 上式より 𝑥1 𝑘 + 1 𝑥2 𝑘 + 1 1 = 0 1 − 𝑝 Τ𝑏 2 𝑥1 𝑘 𝑝 𝑥2 𝑘 + 𝑝 + 𝑏𝑇 − 1 Τ𝑏 2 𝑢(𝑘) 1−𝑝 39

8.4 模型追従適応制御(MRAC) 𝑦 𝑘 + 1 − 𝑦ො 𝑘 + 1 目標： 40 摩擦係数が時間と共に 変わるサーボ機構 → 0 となる入力 u(k) を求める 模型適合則： 𝑢 𝑘 = ① 1 −𝑭𝑇0 ∙ 𝑾0 𝑘 + 𝑩𝑇 ∙ 𝑹(𝑘) 𝑐0 𝑭𝑇0 = 𝑐1 𝑑1 + 1 + 𝑝 𝑾𝑇0 𝑘 = 𝑢 𝑘 𝑦 𝑘 ② −𝑝 ③ 𝑦 𝑘−1 𝑩 ∈ ℝ3×1 , 𝑹 ∈ ℝ3×1 , 𝑭𝟎 ∈ ℝ3×1 , 𝑾𝟎 ∈ ℝ3×1 (8-32) 指令値 逐次での推定アルゴリズム ① 1 ෡ 𝑇0 (𝑘) ∙ 𝑾0 𝑘 + 𝑩𝑇 (𝑘) ∙ 𝑹(𝑘) −𝑭 𝑐0Ƹ (8-35) 新たな入力u(k+1) 𝑢 𝑘 = を計算 ෡ 𝑘+1 =𝑭 ෡ 𝑘 + 𝑮 ∙ 𝑒 ∗ (𝑘 + 1) ∙ 𝑾 𝑘 ③ 𝑭 適応アルゴリズムにより 推定値F(k+1)を計算 (8-36) 入力を加えて 誤差E(z)を算出 𝑇 ② ෡ 𝑘 𝑭−𝑭 𝑾(𝑘) ∗ 𝑒 𝑘+1 = 1 + 𝑾𝑇 (𝑘)𝑮𝑾(𝑘) (8-40) ෡ ∈ ℝ4×1 , 𝑭 ෡𝟎 ∈ ℝ3×1 𝑾 ∈ ℝ4×1 , 𝑮 ∈ ℝ1×1 𝑭

8.4 模型追従適応制御（例） 41 0.5≦ b≦1.5 固定ゲインの 有限整定制御 ⚫ b(k) (粘弾性係数)の変化幅 は，y(k)出力挙動に変化を もたらすようだ ⚫ 固定ゲインの有限整定制御 ではovershootや発振が起 こる --y1(k)-⚫ このシミュレーションでは 外乱の影響がほぼ出ない b(k)：粘弾性係数 scalex40 y1[k]:固定ゲインの有限整定制御 入力u(k)の変化せわしない．こんな ものなのででしょうか？ 書籍List8-2 変更 line330 W(4)→W(3) DGC_no6/p166v1

8.4 模型追従適応制御（例） 42 外乱の影響 固定ゲインの 有限整定制御 ⚫ 外乱をV=300とした 結果を示す．この程度 だとその影響が確認で きる b=1 b=0.5 前頁図はV=50のstep状 の外乱．この程度では大 きな変化が見られない． 入力uの値と比べると納得 𝑢 𝑘 = 1 ෡ 𝑇0 (𝑘) ∙ 𝑾0 𝑘 + 𝑩𝑇 (𝑘) ∙ 𝑹(𝑘) −𝑭 𝑐0Ƹ u(k)が大きな値をとるのはc0_hat の値が小さいため(0.005程度） DGC_no6/p166v1

8.4 模型追従適応制御（例） 書籍の図 page166 “高橋安人,ディジタル制御 ,岩波書店,1987 43

式の導出(MRAC) プラント： 𝐴 𝑧 −1 𝑦 𝑘 = 𝑧 −1 ∙ 𝐶 𝑧 −1 ∙ 𝑢 𝑘 44 (8-23) 𝑦 𝑘 𝐴 𝑧 −1 = 1 − 𝑝𝑧 −1 1 − 𝑧 −1 = 1 − 1 + 𝑝 𝑧 −1 + 𝑝𝑧 −2 𝐶 𝑧 −1 = 𝑐0 + 𝑐1 𝑧 −1 同じ形式とする 固定モデル： 𝐴መ 𝑧 −1 𝑦ො 𝑘 = 𝑧 −1 ∙ B ෡ 𝑧 −1 ∙ 𝑟 𝑘 (8-22)より 𝑌(𝑧) = 𝑧 −1 𝑏0 + 𝑏1 𝑧 −1 𝑅(𝑧) 𝑦ො 𝑘 (8-25) 𝐴መ = 1 𝐵෠ = 𝑏0 + 𝑏1 𝑧 −1 (8-26) プラントとモデル の一致条件： 出力差： 1 + 𝑑1 𝑧 −1 𝑦 𝑘 + 1 − 𝑦ො 𝑘 + 1 =0 d1はなにか？ ちょっと天下り 𝑒 𝑘 = 𝑦 𝑘 − 𝑦ො 𝑘 1 + 𝑑1 𝑧 −1 𝑒 𝑘 + 1 = 𝑒 𝑘 + 1 + 𝑑1 𝑒 𝑘 = 0 𝑒 0 ≠ 0 → 𝑒 1 = 𝑑1 𝑒(0) → 𝑒 2 = 0 なぜ 𝑑1 が必要かは後の話し (8-24) 𝑑1 一般形 𝐷 𝑧 −1 = 1 + 𝑑1 𝑧 −1 + 𝑑2 𝑧 −2 + ⋯ (8-27)

式の導出(MRAC) プラント 8-23は： 𝐴 𝑧 −1 𝑧𝑦 𝑘 = 𝐶 𝑧 −1 ∙ 𝑢 𝑘 これを8-27に加える： 45 −𝐴 𝑧 −1 𝑦 𝑘 + 1 + 𝐶 𝑧 −1 𝑢 𝑘 = 0 𝐶 𝑧 −1 𝑢 𝑘 + 𝐷 𝑧 −1 − 𝐴 𝑧 −1 𝑦 𝑘 + 1 − 𝐷 𝑧 −1 𝑦ො 𝑘 + 1 𝐷 𝑧 −1 − 𝐴 𝑧 −1 = 1 + 𝑑1 𝑧 −1 − 1 − 1 + 𝑝 𝑧 −1 + 𝑝𝑧 −2 = 𝑑1 − 1 + 𝑝 (8-28) (8-29) 𝑧 −1 − 𝑝𝑧 −2 だから モデル8-25は： ෡ 𝑧 −1 𝑟 𝑘 𝐷 𝑧 −1 𝑦ො 𝑘 + 1 = 𝐷 𝑧 −1 𝑧𝑦ො 𝑘 = 𝐷 𝑧 −1 B サーボの場合は： 𝐷 𝑧 −1 𝑦ො 𝑘 + 1 = 1 + 𝑑1 𝑧 −1 𝑦ො 𝑘 + 1 = 1 + 𝑑1 𝑧 −1 𝑧 −1 𝑏0 + 𝑏1 𝑧 −1 𝑟 𝑘 = 𝑏0 𝑟 𝑘 + 𝑏1 + 𝑑1 𝑏0 𝑟 𝑘 − 1 + 𝑑1 𝑏1 𝑟 𝑘 − 2 (8-22)を用いて 𝐴መ = 1 これをベクトル形式でまとめると 𝐷 𝑧 −1 𝑦ො 𝑘 + 1 = 𝑩𝑇 ∙ 𝑹(𝑘) ベクトルで 表現すると ただし 𝑩𝑇 = 𝑏0 𝑏1 + 𝑑1 𝑏0 𝑹𝑇 (𝑘) = 𝑟 𝑘 𝑟 𝑘−1 𝐷 𝑧 −1 = 1 + 𝑑1 𝑧 −1 𝑑1 𝑏1 (8-30) 𝑟 𝑘−2 𝑩 ∈ ℝ3×1 , 𝑹 ∈ ℝ3×1

式の導出(MRAC) 8-28から： 𝐶 𝑧 −1 𝑢 𝑘 + 𝐷 𝑧 −1 − 𝐴 𝑧 −1 46 𝑦 𝑘 + 1 − 𝐷 𝑧 −1 𝑦ො 𝑘 + 1 (8-29) (8-30) Eq.8-28にEq.8-29,30を代入して整理すると 𝑐0 + 𝑐1 𝑧 −1 𝑢 𝑘 + (𝑑1 + 1 + 𝑝) −𝑝 𝑦(𝑘) − 𝑩𝑇 𝑹 𝑘 = 0 𝑦(𝑘 − 1) (8-31) これをベクトル形式でまとめると 模型適合則： 𝑢 𝑘 = 1 −𝑭𝑇0 ∙ 𝑾0 𝑘 + 𝑩𝑇 ∙ 𝑹(𝑘) 𝑐0 𝑭𝑇0 = 𝑐1 𝑑1 + 1 + 𝑝 𝑾𝑇0 𝑘 = 𝑢 𝑘 𝑦 𝑘 −𝑝 𝑦 𝑘−1 𝑭𝟎 ∈ ℝ3×1 , 𝑾𝟎 ∈ ℝ3×1 (8-32)

式の導出（ MRAC時変推定形） 実際に使う模型規範適応制御(MRAC) 47 (8-32)を基に変形 仮定：プラントのパラメータc0が未知．摩擦係数bは時間と共にゆっくり変化する 𝑭T = 𝑐0 𝑭𝑇0 = 𝑐0 𝑐1 𝑑1 + 1 + 𝑝 𝑭 ∈ ℝ4×1 −𝑝 ෡ ∈ ℝ4×1 𝑭 (k) 時変推定形： 逐次での推定 アルゴリズム① ෡ 𝑇 (𝑘) = 𝑐0Ƹ 𝑭 (8-33) ෡ 𝑇0 (𝑘) = 𝑐0Ƹ (𝑘) 𝑐1Ƹ (𝑘) 𝑭 −𝑝(𝑘) Ƹ (8-34) 𝑩 ∈ ℝ3×1 , 𝑹 ∈ ℝ3×1 , (8-35) 𝑑1 + 1 + 𝑝(𝑘) Ƹ 1 ෡ 𝑇0 (𝑘) ∙ 𝑾0 𝑘 + 𝑩𝑇 (𝑘) ∙ 𝑹(𝑘) 𝑢 𝑘 = −𝑭 𝑐0Ƹ 𝑭𝟎 ∈ ℝ3×1 , 𝑾𝟎 ∈ ℝ3×1 推定パラメータを求める適応アルゴリズム 逐次での推定 アルゴリズム③ ෡ 𝑘+1 =𝑭 ෡ 𝑘 + 𝑮 ∙ 𝑒 ∗ (𝑘 + 1) ∙ 𝑾0 𝑘 𝑭 𝑾T 𝑘 = 𝑢 𝑘 𝑾𝑇0 𝑘 = 𝑢 𝑘 ෡ 𝑘+1 𝑒∗ 𝑘 + 1 = 𝑭 − 𝑭 ෡ 𝑘 𝐷 𝑧 −1 𝑒 𝑘 + 1 = 𝑭 − 𝑭 𝑇 𝑮: Gain 𝑢 𝑘−1 𝑦 𝑘 𝑦 𝑘−1 𝑾(𝑘) (8-37) (8-38) 𝑾(𝑘) 𝑇 (8-36) ←証明省く ෡ ∈ ℝ4×1 . 𝑾𝟎 ∈ ℝ3×1 , 𝑾 ∈ ℝ4×1 . 𝑮 ∈ ℝ1×1 𝑭 (8-39)

式の導出（ MRAC時変推定形） 48 Eq. 8-38 右辺に(k+1)項を含む．そこでEq.8-15と似た関係を用いて， 𝑇 ෡ 𝑘 +𝑭 ෡ 𝑘 −𝑭 ෡ 𝑘+1 𝑒∗ 𝑘 + 1 = 𝑭 − 𝑭 𝑾(𝑘) よって 8-36を代入： 逐次での推定 アルゴリズム② 𝑒∗ 𝑘 + 1 ∗ ෡ 𝑘 + 𝑮 ∙ 𝑒 ∗ (𝑘 + 1) ∙ 𝑾0 𝑘 = 𝑭−𝑭 𝑒 𝑘+1 = ෡ 𝑘 𝑭−𝑭 𝑇 𝑇 𝑾(𝑘) ෡ ∈ ℝ4×1 , 𝑾 ∈ ℝ4×1 𝑩 ∈ ℝ3×1 , 𝑭 𝑾(𝑘) 1 + 𝑾𝑇 (𝑘)𝑮𝑾(𝑘) (8-40) サーボ系の例では： 𝑭𝑇 𝑾 𝑘 = y k + 1 + 𝑑1 y k 逐次での推定 アルゴリズム② 𝑒∗ 𝑘 + 1 y k + 1 + 𝑑1 y k − 𝑭𝑇 𝑘 𝑾(𝑘) = 1 + 𝑾𝑇 (𝑘)𝑮𝑾(𝑘) (8-41) ෡ ∈ ℝ4×1 , 𝑾𝟎 ∈ ℝ 3 × 1�, 𝑾 ∈ ℝ 4 × 1�, 𝑮 𝑭 ∈ ℝ4×1 , 𝑭 ∈ ℝ1×1 アルゴリズム 完成

8.4 模型追従適応制御（List8-2） 49 書籍のlist page167 書籍List8-2 変更 line330 W(4)→W(3) 書籍line330の記述だと y(k)→y(k-1) となる たぶん line330 W(3)? 𝑒∗ 𝑘 + 1 y k + 1 + 𝑑1 y k − 𝑭𝑇 𝑘 𝑾(𝑘) = 1 + 𝑾𝑇 (𝑘)𝑮𝑾(𝑘) (8-41) “高橋安人,ディジタル制御 ,岩波書店,1987

8.4 模型追従適応制御（書籍List8-2 変更なし） 0.46≦ b≦1.46 書籍List8-2 変更なし line330 W(4) ｂ：粘弾性係数 50 0.5≦ b≦1.5 修正なしでは：ｂの増加でovershootが発生する（修正済みでは発生しない） DGC_no6/p166v1

8.5 最小2乗同定アルゴリズム 8.1 簡単な適応サーボ系 8.2 模型規範適応制御 8.3 MRASによる同定と適応制御への応用 8.4 模型追従適応制御 8.5 最小2乗同定アルゴリズム 最も基本的なアルゴリズム である最小二乗法の解説 51

8.5 最小2乗同定アルゴリズム 52 最も基本的なアルゴリズムである最小二乗法をとりあげる 同定対象： 𝑎1 , 𝑎2 , ⋯ , 𝑎2 , 𝑏1 , 𝑏2 , ⋯ , 𝑏2 伝達関数の係数を求める Y(z) 𝑏1 𝑧 𝑛−1 + 𝑏2 𝑧 𝑛−2 + ⋯ + 𝑏𝑛−1 𝑧 + 𝑏𝑛 𝑏1 𝑧 −1 + 𝑏2 𝑧 −2 + ⋯ + 𝑏𝑛 𝑧 −𝑛 𝐺 𝑧 = = = U(z) 𝑧 𝑛 + 𝑎1 𝑧 𝑛−1 + 𝑎2 𝑧 𝑛−2 + ⋯ + 𝑎𝑛−1 𝑧 + 𝑎𝑛 1 + 𝑎1 𝑧 −1 + 𝑎2 𝑧 −2 + ⋯ + 𝑎𝑛 𝑧 −𝑛 (8-42) 例としてn=2のとき，次のようになる 𝑦 𝑗 = −𝑎1 𝑦 𝑗 − 1 − 𝑎2 𝑦 𝑗 − 2 + 𝑏1 𝑢 𝑗 − 1 + 𝑏2 𝑢 𝑗 − 2 = 𝒇𝑇 ∙ 𝒙(𝑗) 𝒇𝑇 = −𝑎1 , −𝑎2 , 𝑏1 , 𝑏2 𝒙𝑇 = 𝑦 𝑗 − 1 , 𝑦 𝑗 − 2 , 𝑢(𝑗 − 1), 𝑢(𝑗 − 2) (8-43) 推定値 𝒇෠ の同定： 𝑦 𝑗 − 𝒇෠ 𝑇 ∙ 𝒙 𝑗 = 𝑒(𝑗) (8-44) 評価関数： (8-45) 𝐽 = 𝑤 𝑘−1 𝑒 2 1 + 𝑤 𝑘−2 𝑒 2 2 + ⋯ + 𝑤𝑒 2 𝑘 − 1 + 𝑒 2 𝑘 wが１より小さいときは古いデータへの重みを減らす（過去を忘れる）

8.5 最小2乗同定アルゴリズム（バッチ方式） 53 評価関数Jを最小化するアルゴリズム 𝑦(1) 𝑥 𝑇 (1) 𝑒(1) 𝑇 (2) 𝑦(2) 𝑥 ෠ 𝑒(2) = ∙ 𝒇+ ⋮ ⋮ ⋮ 𝑦(𝑘) 𝑒(𝑘) 𝑥 𝑇 (𝑘) 𝒀(𝑘) 𝑿(𝑘) 𝑬(𝑘) 𝑤 𝑘−1 0 𝑾(𝑘) = ⋮ ⋮ 0 ⋯ ⋱ ⋯ ⋯ ෠ 𝒀 𝑘 = 𝑿(𝑘) ∙ 𝒇𝑬(𝑘) 評価関数： (k)を省略 最小化条件： 最適推定値： ⋯ 𝑤 𝑘−2 𝑤 0 0 ⋮ ⋮ 0 1 対角行列 (8-46) 𝑇 𝐽 𝑘 = 𝑬𝑇 𝑘 𝑾 𝑘 𝑬 𝑘 = 𝒀 𝑘 − 𝑿 𝑘 𝒇෠ 𝑾 𝑘 𝑬(𝑘) 𝒀 𝑘 − 𝑿 𝑘 𝒇෠ (8-48) 𝑇 ෠ 𝒇෠ 𝑻 𝑿𝑻 𝑾𝑿𝒇෠ 𝐽 = 𝒀 − 𝑿𝒇෠ 𝑾𝑬 𝒀 − 𝑿𝒇෠ = 𝒀𝑻 𝑾𝒀 − 𝒇෠ 𝑻 𝑿𝑻 𝑾𝒀 − 𝒀𝑻 𝑾𝑿𝒇+ 𝜕𝐽Τ𝜕𝒇෠ = −2𝒀𝑻 𝑾𝑿 + 2𝒇෠ 𝑻 𝑿𝑻 𝑾𝑿 = 𝟎 𝑿𝑻 𝑾𝒇෠ = 𝑿𝑻 𝑾𝒀 𝒇෠ = 𝑿𝑻 (𝑘)𝑾(𝑘)𝑿(𝑘) −1 𝑿𝑻 𝑘 𝑾 𝑘 𝒀(𝑘) 以上は，逐次推定方式ではなく，データが蓄積されたのち一挙に計算する バッチ方式 (8-49)

8.5 最小2乗同定アルゴリズム（反復方式） 54 オンライン同定に向いた反復方式(recursive method) kまでのデータに(k+1) を付け加える 𝑿 𝑘+1 = 𝑿 𝑘 , 𝑥 𝑇 (𝑘 + 1) 𝒀 𝑘+1 = 𝒀 𝑘 𝑦 (𝑘 + 1) (a) ෠ + 1) = 𝑿𝑻 (𝑘 + 1)𝑾(𝑘 + 1)𝑿(𝑘 + 1) −1 𝑿𝑻 𝑘 + 1 𝑾 𝑘 + 1 𝒀(𝑘 + 1) (8-50) (8-49)から 𝒇(𝑘 H(k)の導入 𝑯−1 𝑘 = 𝑿𝑻 𝑘 𝑾 𝑘 𝑿 𝑘 , 𝑯−1 𝑘 + 1 = 𝑿𝑻 𝑘 + 1 𝑾 𝑘 + 1 𝑿 𝑘 + 1 (b) (8-49)→ ෠ 𝒇(𝑘) = 𝑯(𝑘) 𝑿𝑻 𝑘 𝑾 𝑘 𝒀(𝑘) (c) (8-50) → ෠ + 1) = 𝑯(𝑘 + 1) 𝑿𝑻 𝑘 + 1 𝑾 𝑘 + 1 𝒀(𝑘 + 1) 𝒇(𝑘 (d) (a), (b)から 𝑯−1 𝑘 + 1 = 𝑿𝑻 𝑘 + 1 𝑾 𝑘 + 1 𝑿 𝑘 + 1 = 𝑤𝑿𝑻 𝑘 𝑾 𝑘 𝒀 𝑘 + 𝒙 𝑘 + 1 ∙ 𝒙𝑻 𝑘 + 1 = 𝑤𝑯−1 𝑘 + 𝒙 𝑘 + 1 ∙ 𝒙𝑻 𝑘 + 1 (e)

8.5 最小2乗同定アルゴリズム（反復方式） 55 逆行列則(matrix inversion lemma)の適用 (e)は次の ようになる Hの逆行列を 避ける 𝑯 𝑘 + 1 = 𝑤𝑯−1 𝑘 + 𝒙 𝑘 + 1 ∙ 𝒙𝑻 𝑘 + 1 = −1 1 1 1 1 𝑯 𝑘 − ∙ 𝑯 𝑘 ∙ 𝒙 𝑘 + 1 ∙ 𝒙𝑻 𝑘 + 1 ∙ 𝑯 𝑘 ∙ ∙ 𝑤 𝑤 𝑤 𝐷 1 𝐷 = 1 + 𝒙𝑻 𝑘 + 1 𝑯 𝑘 𝒙 𝑘+1 𝑤 (8-51) ෠ + 1) の最適推定を与える最終形 係数ベクトル式 𝒇(𝑘 係数ベクトル式の 逐次推定式 𝒇෠ 𝑘 + 1 = 𝒇෠ 𝑘 + 𝑮(𝑘 + 1) ∙ 𝑦 𝑘 + 1 − 𝑦ො ° (𝑘 + 1) 1 𝑮 𝑘+1 = 𝑯 𝑘 𝒙 𝑘+1 , 𝑤𝐷 𝑦ො ° 𝑘 + 1 = 𝒙𝑻 𝑘 + 1 ∙ 𝒇෠ 𝑘 これから（８-51）は下記式をもちいて更新される H(k)式の 逐次推定式 𝑯 𝑘+1 = 1 1 𝑯 𝑘 − ∙ 𝑮 𝑘 + 1 𝒙𝑻 𝑘 + 1 ∙ 𝑯 𝑘 𝑤 𝑤 1 = 𝑰 − 𝑮 𝑘 + 1 𝒙𝑻 𝑘 + 1 𝑯 𝑘 𝑤 (8-52)

8.5 最小2乗同定アルゴリズム（例） 56 CASE１ k step= 10 F0= [[ 1.59523973] [-0.7952155 ] [ 0.40067768] [ 0.60046925]] DGC_no6/p174v1

8.5 最小2乗同定アルゴリズム（例） 57 CASE２ k step= 10 F0= [[ 1.59971393] [-0.79969521] [ 0.39993237] [ 0.59991643]] DGC_no6/p174v1

8.5 最小2乗同定アルゴリズム（例） 58 Plant Step response u(k)=0.2 k= 0, y= 0 k= 1, y= 0.08 k= 2, y= 0.328 k= 3, y= 0.6608 k= 4, y= 0.99488 k= 5, y= 1.2632 k= 6, y= 1.4252 k= 7, y= 1.4697 k= 8, y= 1.4114 k= 9, y= 1.2825 k= 10, y= 1.1229 k= 11, y= 0.97059 k= 12, y= 0.85464 k= 13, y= 0.79095 k= 14, y= 0.78182 k= 15, y= 0.81814 k= 16, y= 0.88357 k= 17, y= 0.9592 k= 18, y= 1.0279 k= 19, y= 1.0772 k= 20, y= 1.1013 k= 21, y= 1.1002 k= 22, y= 1.0794 k= 23, y= 1.0468 k= 24, y= 1.0114

8.5 最小2乗同定アルゴリズム（例） 59 書籍図8-8にあわせた表示にした -a1,-a2,-b1,-b2表示とした 入力u(k)は乱数を用いているため，書籍とは収束 過程は同じではない．但し，収束値は等しい p174v1_op1

8.5 最小2乗同定アルゴリズム（例） page174 “高橋安人,ディジタル制御 ,岩波書店,1987 60

おわりに 8章終わりました。予定完了です。 • 書籍の奥付をみると著者が73歳ごろ出版されものです。著者の研究に向けた真摯 な態度と情熱に感銘を受けました。今回、pythonで書き直し、内容をまとめてみ たのは、私も同年代となり少しは著者を手本に心がけたいという思いからでした。 • 「昭和62年12月24日 輪講終了」との記入が本書最後の頁にありました。この 当時は様々な勉強会があったので、本書を誰と行ったかなどの記憶が一切ありま せん。ただ少々難解であったこと、果たして使える技術なのだろうか？と思った ような記憶があります。 • しかし、現在Digital技術が進歩し、演算処理が高速かつ驚くほど安価に行えるよ うになりました。またディジタル系では理論的関係が、ほとんどそのままの形で コンピュータのプログラムになるのは分かり易いです。まさに著者の予測のよう にDigital制御が主力となってきました。 • 最後に本書序論を次頁に記します。40年ほど以前の文章ですが本質をとらえてい るのではないかと思えます。 （Pythonのprogramは下記のGitHubに格納してあります） https://github。com/m4881shimojo/YTakahashi_Digital-control/tree/main 61

• https://github

ディジタル制御 序論 1.1 ディジタル制御工学 制御技術はアナログ方式からディジタル方式へと急速に移りつつある．このため にコントロール・エンジニアは次の課題に直面する．それは従来のように検出端や 操作端のハードウェアに明るいことばかりではなく，マイクロプロセッサ技術とそ の製品の躍進を常につかんでいることが要求される．インターフェイス技術もわき まえなければならない． これらと同程度またはそれ以上に大切なのはソフトウェアを開発する能力である． 適切なディジタル制御アルゴリズムによってチップをフルに使いこなすべきである． 折角ディジタル方式を導入しても，アルゴリズムがアナログ方式のコピーでは， ディジタル制御の真価が発揮できない． 本書の内容はディジタル制御論である．在来の制理論はアナログ系（連続時間 系）を上位におく姿勢であった．理論上はその方が微積分によりすっきりした形に なる．しかしディジタル・ハードウェアとは結びつかない．ディジタル系（離散時 間系）では理論的関係が，ほとんどそのままの形でコンピュータのプログラムに書 かれる．そこで本書では理論と平行してプログラムリストをかかげ，机上のパーソ ナル・コンピュータを使いながら学習が進められるようにする. 引用：“高橋安人，ディジタル制御， p.1 ，岩波書店， 1987 ，第4刷” 62

ふろく：HP-Basic 63 HP-Basicの行列演算は コマンドが独特なので 分かりにくいです。 例えば， MAT F0=F(2:4) しかし，心配無用です。 現在は，chatGPTで pythonに変換してく れます この様なことも可能です 1. HP-basic Listを携帯 で撮影 2. 携帯アプリで自動的 に文字情報に変換 3. chatGPTにpythonへ の変換依頼 今回の書籍Listでは一寸問題 もありますが..

64