ディジタル制御（1987年）その１v2 pythonで記述してみる（2. z変換、3. 線形離散時間形）

2024.1.27 ディジタル制御（1987年）その１v2 pythonで記述してみる（2. z変換、3. 線形離散時間形） https://github.com/m4881shimojo/YTakahashi_Digital-control/tree/main 下 条 誠 電気通信大学名誉教授 https://researchmap.jp/read0072509/ “高橋安人，ディジタル制御 ，岩波書店， 1987 The University of Electro-Communications https://www.docswell.com/user/m_shimojo Department of Mechanical Engineering and Intelligent System

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

目 次 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 1990年ごろ輪講で使った本書を眺めていたら懐かしくなり読み始めました。本書のプロ グラムは行列演算を行えるHP-Basicで記述されています。今回、これをPythonで再記述し てみました。 本書では、基本的なマトリクス演算を利用して、状態空間系から同伴形への変換、PID制 御、線形２乗最適制御、同定および適応制御、ランダム信号系などの高度な内容までをプロ グラムを示しながら解説してあります。ただし、その内容は分かり易いものでなく、初学者 が学ぶには無理があると思います。本書の面白いところは、全てのアルゴリズムが、HPBasi言語と付属の基本行列演算で作られている点です。正準形変換の方法やリカチ行列の解 法などの記述が多々あり大変興味深いものです。これらは現在ではpython-controlを使えば 簡単に解けます。しかし、その解法の基本原理まで踏み込み、意外と簡単なアルゴリズムで 計算できることを実感できるのは楽しいです。 本スライドは単に本人の再学習のためのもので、内容メモとプログラムの再記述だけです。 ディジタル制御の解説ではありません。 また、pythonのプログラムは数多くのミスがあります。書籍の数値とはほぼ同じ値であ ることはチェックしましたが、抜けがあると思います。その点ご注意ください。 （Pythonのprogramは下記のGitHubに格納してあります） https://github.com/m4881shimojo/YTakahashi_Digital-control/tree/main

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

1. 1.1 緒 言 ディジタル制御工学 1.2 ディジタル制御系の構成 1.3 量子化、サンプリングおよびホールド 1.4 例－ロボットアームのデッドビート応答 1.5 プログラム作成入門 4

1.1 ディジタル制御工学 本書の内容はディジタル制御論である．在来の制御理論はアナログ系（連続時間 系）を上位に置く姿勢であった．理論上はその方が微積分によりすっきりした形 になる．しかしディジタル・ハードウェアとは結び付かない．ディジタル系（離 散時間系）では理論的関係が，ほとんどそのままの形でコンピュータのプログラ ムに書かれる．そこで本書では理論と平行してプログラミングリストをかかげ， 机上のパーソナルコンピュータを使いながら学習が進められるようにする． （原文そのまま） 引用：“高橋安人，ディジタル制御， p.1 ，岩波書店， 1987 ，第4刷” 5

1.2 ディジタル制御系の構成 v(t) u(t) y(t) y(k) r(k) ① ディジタル演算素子は目標値と制御量の差を消し去るような制御信号を算定する e(k)=r(k)-y(k) ② フィードバック制御系の性能は制御法則に大きく支配される。アナログ制御では電子回路や 空気圧機素などを使って実現する ③ ただし、このハードウェアに拘束され単純な制御法則（PIDなど）しか実現できなかった ④ ディジタル制御では、例えばディジタル演算素子にプラントの数学模型を内蔵させ、実際に は測れないプラントの内部状態を推測し、フィードバックすることができる ⑤ これにより従来よりきめの細かい制御が実現できる 6

１.4 ロボットアームのデッドビート応答 1次系 自由応答 𝑑𝑥 Τ𝑑𝑡 = 𝑎𝑥 𝑥 𝑡 = 𝑥0 𝑒 𝑎𝑡 ( 𝑥0 , 𝑥1 , ⋯ 𝑥𝑘 ⋯) 等比数列 (𝑥𝑘 = 𝑝𝑘 𝑥0 , 𝑝 = 𝑒 𝑎𝜏 ) τ：sampling time リスト１－１ ロボットアームの有限時間整定応答 単位パルス入力uに対する応答ｇと仮定する。 𝑔 𝑘 = 𝑝𝑘−1 , (𝑝 = −0.5など） (1-3) 応答yは次のようになる。 𝑘 𝑦 𝑘 = ෍𝑢 𝑗 ∙ 𝑔 𝑘 − 𝑗 (1-5) 𝑗=0 問題：ロボットアームの位置y(k)を一 定y(k)=1にとどめるため、制御入力 u(k)をどのように設定すればよいか？ そのプログラムと、制御結果の図を次 頁に示す 図 ｐによるg(k)の推移 DGC_no1/ p8v1 7

１.4 ロボットアームのデッドビート応答 デッドビート応答でのu(k)の決定方法 𝑘 𝑦 𝑘 = ෍𝑢 𝑗 ∙ 𝑔 𝑘 − 𝑗 要求：y(k)を0→1とし、そのまま”1”に留める 最初の入力u(0)： 𝑦 0 = 𝑢(0) ∙ 𝑔(0) 𝑗=0 ex.（P=-0.5とする） g(0)=0だからy(0)=0 次の入力 u(1)： 𝑦 1 = 𝑢 0 ∙ 𝑔 1 + 𝑢(1) ∙ 𝑔(0) g(0)=0, g(1)=1でy(1)となるためには → 𝒖 𝟎 = 𝟏 次の入力 u(2)： 𝑦 2 = 𝑢 0 ∙ 𝑔 2 + 𝑢(1) ∙ 𝑔(1) +𝑢(2) ∙ 𝑔(0) g(0)=0, g(1)=1, g(2)=-0.5でy(2)=1となるためには、 𝑦 2 = 1 ∙ (−0.5) + 𝑢(1) ∙ 1 +𝑢 2 ∙ 0 = 1 -1 ∙ 5 + 𝑢(1) = 1 この計算を反復すると、k>0に対してu(k)=1.5 → 𝒖 𝟏 = 𝟏. 𝟓 とすればよい 8

１.4 ロボットアームのデッドビート応答 9 図1-3 ロボットアーム の有限時間整定応答 (P=-0.5) DGC_no1/ page10

１.4 ロボットアームのデッドビート応答 10 図1-3 ロボットアーム の有限時間整定応答 (P=0.4) DGC_no1/ page10

2.1 ｚ変換と反復計算 2.1 ｚ領域における時系列 2.2 入力と出力 2.3 極とゼロ 2.4 調和振動 2.5 リカーシブ・アルゴリズム 11

2.1 ｚ領域における時系列 12 ラプラス変換の離散時間バージョンがｚ変換である 𝑧 −1 は１サンプリング間隔（周期T）にあたる時間遅れを表す演算子である。 −𝑠𝑇 ラプラス変換のパラメータｓを使うと、この時間遅れ（むだ時間）は 𝑒 となる （ｚはＴだけの時間を進めることを意味する） 一般に、k=0,1,2..での信号値がf(0), f(1), f(2),..で表される時系列 のz変換は次のように定義される 時系列信号： 𝑍 = 𝑒 𝑎𝑇 𝑓 0 ,𝑓 1 ,𝑓 2 ,𝑓 3 ,⋯ ∞ z変換： 𝒵 𝑓(𝑘) = 𝑓 0 + 𝑓 1 ∙ 𝑧 −1 + 𝑓 2 ∙ 𝑧 −2 ⋯ = ෍ 𝑓 𝑘 ∙ 𝑧 −𝑘 𝑘=0 x 𝑓 𝑡 x 連続信号 𝑓 1 サンプル信号 𝑓 2 𝑓 3 𝑓 0 x0 𝑓 4 T T 2T 3T 4T t t sampler T 2T 3T 4T

2.1 ｚ領域における時系列 例） 𝒵 𝑓(𝑘) = 𝑓 0 + 𝑓 1 ∙ 𝑧 −1 + 𝑓 2 ∙ 𝑧 −2 ⋯ F(z)=1/(z-p)は次のように展開できる 𝐹 𝑧 = 等比時系列 のz変換 1 = 1 ∙ 𝑧 −1 + 𝑝𝑧 −2 + 𝑝 2 𝑧 −3 + 𝑝3 𝑧 −4 ⋯ 𝑧−𝑝 f(0) p=1に とると 𝐹 𝑧 = f(1) 13 f(2) f(3) 1 , p , p 2, p 3… . これは 等比系列での時系列信号 1 = 𝑧 −1 + 𝑧 −2 + 𝑧 −2 + 𝑧 −3 ⋯ 𝑧−1 k=1からの単位 ステップを表す (下図左) 信号値がf(0)=0, f(1)=1, f(2)=1,..で表される時系列 単位 ステップ 𝐹 𝑧 = 𝑧 = 1 + 𝑧 −1 + 𝑧 −2 + 𝑧 −2 + 𝑧 −3 ⋯ 𝑧−1 k=0からの単位 ステップを表す (下図右) k=0からの単位ステップは、zを掛けて1刻み進める zを掛ける p等比系列 とは？ t T 2T 3T 4T t T 2T 3T 4T

2.1 ｚ領域における時系列 14 ｐの等比系列は、1次系の自由応答を表す x x0 x0 連続信号 x x1 サンプル信号 x2 T x3 x4 t t sampler 1次系の自由応答： 𝑑𝑥(𝑡) = 𝑎𝑥(𝑡) 𝑑𝑡 T 2T 3T 4T 𝑥 𝑡 = 𝑒 𝑎𝑡 𝑥0 𝑎𝑇 サンプル信号： 𝑥1 = 𝑒 𝑥0 , 𝑥2 = 𝑒 𝑎(2𝑇) 𝑥0 , 𝑥3 = 𝑒 𝑎(3𝑇) 𝑥0 , 𝑥4 = 𝑒 𝑎(4𝑇) 𝑥0 , ⋯ 等比級数となる 𝑥1 = 𝑝𝑥0 , 𝑥2 = 𝑝2 𝑥0 , 𝑥3 = 𝑝3 𝑥0 , ⋯ 𝑥𝑘 = 𝑝𝑘 𝑥0 , ⋯ 1次系自由応答のサンプル値は等比級数となる ｚ変換形は 𝐹 𝑧 = 𝑝 = 𝑒 𝑎𝑇 1, 𝑝, 𝑝 2 , 𝑝3 , ⋯ 1 = 1 ∙ 𝑧 −1 + 𝑝𝑧 −2 + 𝑝2 𝑧 −3 + 𝑝3 𝑧 −4 ⋯ 𝑧−𝑝

2.2 入力と出力 2.1 ｚ領域における時系列 2.2 入力と出力 2.3 極とゼロ 2.4 調和振動 2.5 リカーシブ・アルゴリズム 15

2.2 入力と出力 16 状態式と出力式は、次の形へ一般化される 差分式 Z変換式 状態式： 𝑥 𝑘 + 1 = 𝑝𝑥 𝑘 + 𝑞𝑢(𝑘) zX z − 𝑧𝑥 0 = 𝑝𝑋 𝑧 + 𝑞𝑈(𝑧) 出力式： 𝑦 𝑘 = 𝑐𝑥 𝑘 + 𝑑𝑢(𝑘) Y 𝑧 = 𝑐𝑋 𝑧 + 𝑑𝑈(𝑧) 𝒵 𝑥(𝑘) = 𝑋 𝑧 , 𝒵 𝑢(𝑘) = 𝑈 𝑧 , 𝒵 𝑦(𝑘) = 𝑌 𝑧 解説： 状態式は、kから(k+1)へ時間変化を表す。このz変換式は、 初期値x(0)を用い、次のようになる 変換式の定義から、 𝒵 𝑥(𝑘) = 𝑥 0 + 𝑥 1 ∙ 𝑧 −1 + 𝑥 2 ∙ 𝑧 −2 ⋯ 定義から 状態変数x(k+1)は、 𝒵 𝑥(𝑘 + 1) = 𝑥 1 + 𝑥 2 ∙ 𝑧 −1 + 𝑥 3 ∙ 𝑧 −2 ⋯ = 𝑧 𝑥 0 + 𝑥 1 ∙ 𝑧 −1 + 𝑥 2 ∙ 𝑧 −2 ⋯ − 𝑧𝑥 0 よって、 𝒵 𝑥(𝑘 + 1) = zX z − 𝑧𝑥 0 X(z)

2.1 入力と出力 17 応答を求めてみる 自由応答： 𝑥 0 ≠ 0, 𝑈 𝑧 = 0 → 𝑧 𝑋 𝑧 = 𝑥(0) 𝑧−𝑝 𝑐𝑧 𝑌 𝑧 = 𝑥(0) 𝑧−𝑝 強制応答： 𝑥 0 = 0, 𝑈 𝑧 ≠ 0 → 𝑞 𝑋 𝑧 = 𝑈(𝑧) 𝑧−𝑝 𝑌 𝑧 = 𝑑𝑧 + 𝑐𝑞 − 𝑝𝑑 𝑈(𝑧) 𝑧−𝑝 自由応答： 𝑧 𝑋 𝑧 = 𝑥(0) = 1 + 𝑝𝑧 −1 + 𝑝 2 𝑧 −2 + 𝑝3 𝑧 −3 ⋯ 𝑥(0) 𝑧−𝑝 x(0), x(1), x(2), x(3) 強制応答： 𝑌 𝑧 = 𝐺(𝑧) ∙ 𝑈(𝑧) 𝐺 𝑧 = 𝑑𝑧 + 𝑐𝑞 − 𝑝𝑑 𝑧−𝑝 入力が単位ステップ U 𝑧 = 𝑧 𝑧−1 𝑥 𝑘 = 𝑝𝑘 𝑥(0) 𝑦 𝑘 = 𝑐 𝑝𝑘 𝑥(0)

2.3 極とゼロ 2.1 ｚ領域における時系列 2.2 入力と出力 2.3 極とゼロ 2.4 調和振動 2.5 リカーシブ・アルゴリズム 18

2.3 極とゼロ 例）2次の伝達関数 z1：ゼロ 𝑧 − 𝑧1 G 𝑧 =𝐾 𝑧 − 𝑝1 𝑧 − 𝑝2 部分分数展開： G 𝑧 = 𝐶1 𝐶2 + 𝑧 − 𝑝1 𝑧 − 𝑝2 19 p1, p2：極、または固有値 (𝑝1 ≠ 𝑝2 ) −1 −2 単位パルス入力： 𝑈 𝑧 = 𝒵 𝑢(𝑘) = 𝑢 0 + 𝑢 1 ∙ 𝑧 + 𝑢 2 ∙ 𝑧 ⋯= 1 1 出力： 𝑌 𝑧 = G 𝑧 𝑈(𝑧) = この時系列出力の和： 0 0 𝐶1 𝐶2 + 𝑧 − 𝑝1 𝑧 − 𝑝2 0, 𝐶1 , 𝐶1 𝑝1 , 𝐶1 𝑝1 2 , ⋯ 0, 𝐶2 , 𝐶2 𝑝2 , 𝐶2 𝑝2 2 , ⋯ 極のモード(mode) 重極のとき： (𝑝1 = 𝑝2 = 𝑝) Y 𝑧 = 𝐶1 𝐶2 + 𝑧−𝑝 𝑧−𝑝 2 3 この時系列出力の和： 0, 𝐶1 , 𝐶1 𝑝 , 𝐶1 𝑝 , 𝐶1 𝑝 , ⋯ 2 0, 0, 𝐶2 , 2𝐶2 𝑝 , 3𝐶2 𝑝2 , 4𝐶2 𝑝3, ⋯ 極のモード(mode)

2.4 調和振動 2.1 ｚ領域における時系列 2.2 入力と出力 2.3 極とゼロ 2.4 調和振動 2.5 リカーシブ・アルゴリズム ⚫ ⚫ ⚫ ⚫ ⚫ リスト2-1 リスト2-2 リスト2-4 リスト2-5 リスト2-6 等比数列 等差数列 平均値の算定 重みを含む平均値の算定 ５点の移動平均 以上のリストは割愛しました 20

2.4 調和振動 21 2次系の状態式と出力式 𝑥1 𝑘 + 1 𝑥2 𝑘 + 1 状態式： 出力式： 𝑦(𝑘) = 𝑐1 𝑝11 = 𝑝 21 𝑐2 𝑝12 𝑝22 𝑥1 𝑘 𝑥2 𝑘 𝑥1 𝑘 𝑥2 𝑘 𝑞1 + 𝑞 u(k) 2 (2-21) + 𝑑u(k) 伝達関数が部分分数展開可能ならば、pijの2×2行列が対角化される。 しかし、極が重複すると対角化ができない。 ⚫ 対角化例 ⚫ 極が重複化例 1 𝐺 𝑧 = 𝑧−𝑝 2 𝑥1 𝑘 + 1 = 𝑝1 𝑥1 𝑘 + 𝑐1 𝑢 𝑘 𝑥2 𝑘 + 1 = 𝑝2 𝑥2 𝑘 + 𝑐2 𝑢 𝑘 𝑦 𝑘 = 𝑥1 𝑘 + 𝑥2 𝑘 (2-22) 𝑥1 𝑘 + 1 = 𝑝 𝑥1 𝑘 + 𝑥2 𝑘 𝑥2 𝑘 + 1 = 𝑝𝑥2 𝑘 + 𝑢 𝑘 𝑦 𝑘 = 𝑥1 𝑘 (2-23)

2.4 調和振動 22 固有値が共役複素数であるシステムのモードは振動となる 𝑝1 = 𝛼 + 𝑗𝛽, 𝑝2 = 𝛼 − 𝑗𝛽 𝑗= −1, 𝛼, 𝛽: 𝑟𝑒𝑎𝑙 𝑛𝑢𝑚𝑏𝑒𝑟 (2-24) 正弦波のz変換: 𝑓 𝑘 = sin 𝑘𝜔𝑇 = 𝑒 𝑗𝑘𝜔𝑇 − 𝑒 −𝑗𝑘𝜔𝑇 Τ 2𝑗 𝐹 𝑧 = 𝑧 −1 𝑒 𝑗𝑘𝜔𝑇 + 𝑧 −2 𝑒 2𝑗𝑘𝜔𝑇 + ⋯ ൗ 2𝑗 − 𝑧 −1 𝑒 −𝑗𝑘𝜔𝑇 + 𝑧 −2 𝑒 −2𝑗𝑘𝜔𝑇 + ⋯ Τ 2𝑗 1 𝑒 𝑗𝑘𝜔𝑇 𝑒 −𝑗𝑘𝜔𝑇 1 𝑧 𝑒 𝑗𝑘𝜔𝑇 − 𝑒 −𝑗𝑘𝜔𝑇 = − = 2𝑗 𝑧 − 𝑒 𝑗𝑘𝜔𝑇 𝑧 − 𝑒 −𝑗𝑘𝜔𝑇 2𝑗 𝑧 2 − 𝑧 𝑒 𝑗𝑘𝜔𝑇 + 𝑒 −𝑗𝑘𝜔𝑇 + 1 𝑧 sin 𝜔𝑇 = 2 (2-27) 𝑧 − 2𝑧 cos 𝜔𝑇 + 1 余弦波のz変換: 𝑓 𝑘 = cos 𝑘𝜔𝑇 = 𝑒 𝑗𝑘𝜔𝑇 + 𝑒 −𝑗𝑘𝜔𝑇 Τ2 𝑧 2 − 𝑧 cos 𝜔𝑇 𝐹 𝑧 = 2 𝑧 − 2𝑧 cos 𝜔𝑇 + 1 (2-28) 根 𝑧 2 − 2𝑧 cos 𝜔𝑇 + 1 = 0 𝑝1 = 𝛼 + 𝑗𝛽, 𝑝2 = 𝛼 − 𝑗𝛽 の形をとる

調和振動 2.4 23 例）発振器モデルを作ってみる 𝛼 = 𝑟 ∙ cos 𝜔𝑇 , 𝛽 = 𝑟 ∙ 𝑠𝑖𝑛 𝜔𝑇 とする 𝑝1 , 𝑝2 = 𝛼 ± 𝑗𝛽 𝑧 2 − 2𝑧 cos 𝜔𝑇 + 1 = 0 𝑥1 𝑘 + 1 𝑥2 𝑘 + 1 = 𝛼 −𝛽 𝑦(𝑘) = 𝑐1 𝑧 2 − 2𝑟𝑧 cos 𝜔𝑇 + 𝑟 2 = 0 𝛽 𝛼 𝑥1 𝑘 𝑥2 𝑘 𝑐2 𝑥1 𝑘 𝑥2 𝑘 𝑞1 + 𝑞 u(k) 2 + 𝑑u(k) DGC_no1/ page24 r=0.95 1周期が12点の 時系列を発生 𝑟 = 1.0 , 𝜔𝑇 = 2𝜋Τ12 𝑟 = 0.95, 𝜔𝑇 = 2𝜋Τ12

調和振動 2.4 24 状態変数x1,x2を以下のようにする 𝑥1 𝑡 = sin 𝜔𝑡 𝛼 = cos 𝜔𝑇 , 𝑥2 𝑡 = cos 𝜔𝑡 𝛽 = sin(𝜔𝑇) そのT秒後（１サンプル時間後）の値は次のようになる 𝑥1 𝑡 + 𝑇 = sin 𝜔 𝑡 + 𝑇 = sin 𝜔𝑡 cos 𝜔𝑇 + cos 𝜔𝑡 cos(𝜔𝑇) 𝑥2 𝑡 + 𝑇 = 𝑐𝑜𝑠 𝜔 𝑡 + 𝑇 = 𝑐𝑜𝑠 𝜔𝑡 cos 𝜔𝑇 − 𝑠𝑖𝑛 𝜔𝑡 𝑠𝑖𝑛(𝜔𝑇) これを状態空間モデルで表現すると次のようになる 𝑥1 𝑡 + 𝑇 𝑥2 𝑡 + 𝑇 = 𝛼 −𝛽 𝛽 𝛼 𝑥1 𝑡 𝑥2 𝑡 𝑞1 + 𝑞 u(t) 2 𝑦 𝑡 = 𝑐1 𝑐2 𝑥1 𝑡 𝑥2 𝑡 例えば、周期をP=2π/ωとし、１周期を12点の時系列で表すとすると、 𝑃 = 2𝜋Τ𝜔 = 12𝑇 ⇒ 𝜔𝑇 = 2𝜋Τ12 となる + 𝑑𝑢(𝑘)

３. 3.1 状態空間系 3.2 伝達関数 3.3 同伴形（正準形） 3.4 応答の算定 3.5 安定に関する定理 線形離散時間系 25

3.1 状態空間系 状態空間系 状態式： 𝒙 𝑘 + 1 = 𝑷𝒙 𝑘 + 𝑸𝒖(𝑘) 出力式： 𝒚 𝑘 = 𝑪𝒙 𝑘 + 𝑫𝒖(𝑘) z変換すると 26 状態式： 𝑧𝑿 𝑧 − 𝑧𝒙(0) = 𝑷𝑋 𝑧 + 𝑸𝑼(𝑧) 出力式： 𝒀 𝑧 = 𝑪𝑿 𝑧 + 𝑫𝑼(𝑧) 自由応答： 𝑿 𝑧 = 𝑧𝑰 − 𝑷 −1 𝑧𝒙 0 𝒙 0 ≠ 0, 𝑼 𝑧 = 0 強制応答： 𝑿 𝑧 = 𝑧𝑰 − 𝑷 −1 𝑸𝑼 𝑧 𝒙 0 = 0, 𝑼 𝑧 ≠ 0 (3-1) (3-3) (3-4) 伝達関数： 𝒀 𝑧 = 𝑮 𝑧 𝑼 𝑧 𝑮 𝑧 = 𝑪 𝑧𝑰 − 𝑷 特性式： −1 𝑸+D (3-5) 𝐹 𝑧 = det 𝑧𝐼 − 𝑃 (3-7) 𝐹(𝑧) = 𝑧 𝑛 + 𝑎1 𝑧 𝑛−1 + 𝑎2 𝑧 𝑛−2 + ⋯ + 𝑎𝑛−1 𝑧 + 𝑎𝑛 (3-8) F(z)=0がシステムの固有値（極）である

3.1 状態空間系（図3-1） 例）4次系の自由応答式 「（2-23）の複素極が重複したシステム」 𝛼 𝛽 −𝛽 𝛼 𝑷= 0 1 0 0 0 0 0 0 𝛼 𝛽 −𝛽 𝛼 𝑪= 0 0 1 0 𝛼 = 𝑟𝑐𝑜𝑠 20° , 𝛽 = 𝑟𝑠𝑖𝑛 20° 0 𝒙(𝟎) = 1 0 ,𝑟 = 1 0 発振器モデル ｒ＝１ 式（2-23）の複素極が重複し たシステムである。z面の単位 円上(r=1)で複素極が重なると、 正弦波のk倍がモードになるた め、時間すなわち刻みkととも に自由応答が成長する。 DGC_no1/ page33 ｒ＝0.98 27 page33

3.2 伝達関数 3.1 状態空間系 3.2 伝達関数 3.3 同伴形（正準形） 3.4 応答の算定 3.5 安定に関する定理 28

3.2 伝達関数 29 伝達関数 B(z) 𝑏1 𝑧 𝑛−1 + 𝑏2 𝑧 𝑛−2 + ⋯ + 𝑏𝑛−1 𝑧 + 𝑏𝑛 𝐺 𝑧 = = A(z) 𝑧 𝑛 + 𝑎1 𝑧 𝑛−1 + 𝑎2 𝑧 𝑛−2 + ⋯ + 𝑎𝑛−1 𝑧 + 𝑎𝑛 𝐺 𝑧 −1 B(𝑧 −1 ) 𝑏1 𝑧 −1 + 𝑏2 𝑧 −2 + ⋯ + 𝑏𝑛−1 𝑧 −(𝑛−1) + 𝑏𝑛 𝑧 −𝑛 = = A(𝑧 −1 ) 1 + 𝑎1 𝑧 −1 + 𝑎2 𝑧 −2 + ⋯ + 𝑎𝑛−1 𝑧 −(𝑛−1) + 𝑎𝑛 𝑧 −𝑛 𝑌(𝑍) = 𝐺(𝑧) 𝑈(𝑘) (3-9) (3-10) A(𝑧 −1 )𝑌(𝑍) = B 𝑧 −1 𝑈(𝑘) z-1は１刻みの遅れを表すから時間領域において以下の式が成り立つ 𝑦(𝑘) +𝑎1 𝑦 𝑘 − 1 + ⋯ + 𝑎𝑛 𝑦 𝑘 − 𝑛 = 𝑏1 𝑢 𝑘 − 1 + ⋯ + 𝑏𝑛 𝑢 𝑘 − 𝑛 𝑦 𝑘 = − 𝑎1 𝑦 𝑘 − 1 + ⋯ + 𝑎𝑛 𝑦 𝑘 − 𝑛 + 𝑏1 𝑢 𝑘 − 1 + ⋯ + 𝑏𝑛 𝑢 𝑘 − 𝑛 (3-11)

3.2 伝達関数 30 システムの動的挙動は、極ばかりでなくゼロによっても大きく左右される page35 例) ステップ応答に及ぼすゼロの影響 G 𝑧 = (type 0) 𝑏1 𝑧 + 𝑏2 , 𝑧 − 0.4 𝑧 − 0.8 z1=0.836 𝑏2 = − 𝑧1 𝑏1 z1=0.6 単位ステップ応答の定常終端状態 𝑏1 + 𝑏2 + ⋯ + 𝑏𝑛 𝑦𝑓 = = 𝐺(1) 1 + 𝑎1 + 𝑎2 + ⋯ + 𝑎𝑛 dc (定常)ゲイン b1=0 z1=1.6 (3-12) dcゲインが1であるようにb1,b2を選んだ STEP応答がk→∞で有限の終端値 G(z)：ゼロ型(type 0)と称する DGC_no1/ page35

3.3 同伴形（正準形） 3.1 状態空間系 本書では同伴形と記述していますが、一般的に 3.2 伝達関数 3.3 同伴形（正準形） 3.4 応答の算定 3.5 安定に関する定理 可制御正準形（Controllable Canonical Form）とは、状態空間表現 を座標変換によって正準形に変換 することで、制御系の解析や設計を 容易にする手法です 可制御正準形と記述されています 31

3.3 同伴形（正準形） 0 0 0 𝑷𝒄 = ⋮ 0 −𝑎𝑛 1 0 0 ⋮ 0 −𝑎𝑛−1 0 1 0 ⋮ 0 ⋯ … 0 ⋯ 0 ⋯ 0 ⋮ ⋮ ⋯ 1 ⋯ −𝑎1 0 0 0 𝒒𝒄 = ⋮ 0 1 32 可制御正準形 𝑪𝒄 = 𝑏𝑛 𝑏𝑛−1 𝑏𝑛−2 ⋯ 𝑏1 伝達関数の分母、分子の係数がPcとCcに明白に含まれる 同伴形は状態空間におけるP, q, cから特性式 F(z) や伝達関数 G(z)=B(z)/A(z), F(z)=A(z) を算定するのに応用される。 𝑝11 ⋯ ⋯ 𝑝1𝑛 ⋯ ⋯ ⋯ ⋯ 𝑷= ⋮ ⋮ ⋮ ⋮ ⋯ ⋯ ⋯ 𝑝𝑛𝑛 0 0 𝑷𝒄 = ⋮ 0 −𝑎𝑛 変換のアルゴリズムを次に記す 1 0 ⋮ 0 −𝑎𝑛−1 … ⋯ ⋮ ⋯ ⋯ 0 0 ⋮ 1 −𝑎1

3.3 同伴形（正準形） （変換方法） 33 可制御同伴形を得るためには変換行列T(n×n)を用いて状態ベクトルを次のように変換する 𝒙𝒄 = 𝑻𝒙 (3-17) 𝒙 𝑘 + 1 = 𝑷𝒙 𝑘 + 𝒒𝑢(𝑘) 𝒙= 𝑻−1 𝒙𝒄 𝑦 𝑘 = 𝑪𝒙 𝑘 𝒙𝒄 𝑘 + 1 = 𝑻𝑷𝑻−1 𝑷𝒙𝒄 𝑘 + 𝑻𝒒 𝑢(𝑘) 𝑦 𝑘 = 𝑪𝑻−1 𝒙𝒄 𝑘 これから、 𝑷𝒄 = 𝑻𝑷𝑻−1 , 𝒒𝒄 = 𝑻𝒒 , 𝑪𝒄 = 𝑪𝑻−1 (3-18) 変換行列Tが 求まればOK 変換行列Tの求め方 Step1:可制御行列とするE(n×n)を作る 𝑬= 𝒒 Step2:システムが可制御の時、E-1の下 端の行（横）ベクトルdを求める 𝒅= 𝟎 𝟎 ⋯ Step3:所要の変換行列Tを、Pとdを用いて 次の形に作成する 𝒅 𝒅𝑷 𝑻 = 𝒅𝑷2 ⋮ 𝒅𝑷𝑛−1 𝑷𝒒 𝑷2 𝒒 ⋯ 1 𝑬−1 𝑷𝑛−1 𝒒 (3-19) (3-20) (3-21) DGC_no1/page40v1

3.3 同伴形（正準形）（プログラム実行結果） 同伴形（正準形）への変換プログラム実行結果 -------------P, Q, C------------------P : [[ 0.921 0.335 0. 0. ] [-0.335 0.921 0. 0. ] [ 0. 1. 0.921 0.335] [ 0. 0. -0.335 0.921]] q : [[0.] [1.] [0.] [0.]] 図3-1の4次系 c : [0. 0. 1. 0.] DGC_no1/page40v1 可制御正準形 ---------Pc,Qc,Cc-------------------An......A1: [ 0.92249 -3.53836 5.3139 -3.684 ] ------------------------------------Bn......B1: [ 0.84824 -1.842 1. 0. ] ------------------------------------Pc : [[ 0. 1. 0. 0. ] [ 0. 0. 1. 0. ] [ 0. 0. 0. 1. ] [-0.92249 3.53836 -5.3139 3.684 ]] qc : [[ 0.] [ 0.] [-0.] [ 1.]] 34 同伴形になる

3.3 同伴形（正準形） 35 Python-Controlを使って検証を行った Python-Controlの出力結果 B(z) 𝐺 𝑧 = A(z) 8.882e-16 s^3 + s^2 - 1.842 s + 0.8481 ---------------------------------------------s^4 - 3.684 s^3 + 5.313 s^2 - 3.538 s + 0.9224 modedeg=20.;r=0.98 a=np.cos(2*np.pi/360*modedeg)*r b=np.sin(2*np.pi/360*modedeg)*r P=np.array([[a,b,0.0,0.0], [-b,a,0.0,0.0], [0.0,1.0,a,b], [0.0,0.0,-b,a]]) 書籍方法での出力結果 Q=np.array([[0], [1.], A: (an,an-1...a2,a1:) [ 0.92249, -3.53836, 5.3139, -3.684 ] ------------------------------------B: (bn,bn-1…b2,b1☺ [ 0.84824, -1.842, 1. , 0. ] [0]]) C=np.array([0,0,1.,0]) # from control import * sys1ss = ss(P, Q, C, 0) sys1tf = ss2tf(sys1ss) print(sys1ss) DGC_no1/page40v1 print(sys1tf) Python-Control [0],

3.3 同伴形（正準形）（プログラム実行結果）36 A:[ 0.92249, -3.53836, 5.3139, -3.684 ] B:[ 0.84824, -1.842, 1., 0. ] 可制御 可観測

3.4 応答の算定 3.1 状態空間系 3.2 伝達関数 3.3 同伴形（正準形） 3.4 応答の算定 3.5 安定に関する定理 37

3.4 応答の算定 38 状態式： 𝒙 𝑘 + 1 = 𝑷𝒙 𝑘 + 𝑸𝒖(𝑘) 出力式： 𝒚 𝑘 = 𝑪𝒙 𝑘 例） 𝛼1 𝛽1 −𝛽1 𝛼1 𝑷= 0 0 0 0 0 0 1 0 𝛼2 𝛽2 −𝛽2 𝛼2 𝑫 = 𝟎 とした 0 𝑞 𝒒= 1 0 𝑞2 𝑪 = 𝑐1 0 𝑐2 0 (3-28) 発振器モデル 固有値α2±jβ2へ別の固有値α1±jβ1の ものを直列につないだ2重振動系である 状態式と出力式の信号線図 各サブシステム の特性式 𝐹1 𝑧 = 𝑧 2 − 2𝛼1 𝑧 + 𝑅1 全体の伝達関数 𝑌(𝑧) 𝛽1 𝑐1 𝑞1 𝐹2 𝑧 + 𝛽2 𝑐2 𝑞2 𝐹1 𝑧 + 𝛽1 𝛽2 𝑐1 𝑞2 𝐺 𝑧 = = 𝑈(𝑧) 𝐹1 𝑧 𝐹2 𝑧 𝑅1 = 𝛼12 + 𝛽12 𝐹2 𝑧 = 𝑧 2 − 2𝛼2 𝑧 + 𝑅2 𝑅2 = 𝛼22 + 𝛽22 (3-29) (3-30)

3.4 応答の算定（可制御と可観測） 場合 入力側 出力側 (a) 𝑞1 =0, 𝑞2 =1 𝑐1 =1, 𝑐2 =0 (b) 𝑞1 =1, 𝑞2 =0 𝑐1 =1, 𝑐2 =0 ←非可制御 (c) 𝑞1 =0, 𝑞2 =1 𝑐1 =0, 𝑐2 =1 ←非可観測 q2=0で切断 𝐹1 𝑧 c1=0で切断 𝐹2 𝑧 サブシステム１ (b)の場合、サブシステム１の状態 変数が入力u(k)の支配を受けない サブシステム２ (c)の場合、サブシステム２の状態 変数が出力y(k)へ関与しない 状態式と出力式の信号線図 （例：2重振動系） 39

3.4 応答の算定（可制御と可観測） 入 力 側 (a)可制御：入力は、q1=0でもq2≠0である限り、サブシステム1からサブシステ ム2へ伝わる。このようにシステムの全ての状態変数x1～x4が入力の 支配を受ける。すなわち可制御である。 (b)非可制御：入力は、q2=0のため、x3, x4左への入力が完全に遮断され、非可 制御になる 出 力 側 (b)可観測：c1≠0である限り、 c2=0でも全ての状態変数x1～x4がの動きが出力に 関係する。すなわち可観測である。 (c)非可観測： c1=0のため、x1, x2の挙動が出力に反映しなくなり、非可観測に おちいる。 (a) 𝐺 𝑧 = 𝛽1 𝛽2 Τ 𝐹1 𝑧 ∙ 𝐹2 𝑧 (b) 𝐺 𝑧 = 𝛽1 𝑐1 𝑞1 𝐹2 𝑧 Τ 𝐹1 𝑧 ∙ 𝐹2 𝑧 =𝛽1 𝑐1 𝑞1 Τ𝐹1 𝑧 (c) 𝐺 𝑧 = 𝛽2 𝑐2 𝑞2 𝐹1 𝑧 Τ 𝐹1 𝑧 ∙ 𝐹2 𝑧 = 𝛽2 𝑐2 𝑞2 Τ𝐹2 𝑧 (b)では非可制御のモードのF2(z)が、 (c)では非可観測のモードのF1(z)が、 何れも入出力関係G(z)から脱落する 極ゼロ相殺 40

3.4 応答の算定（プログラム実行結果） 41 # ケース分けのパラメータ設定 if l==0: q1=0.;q2=1.;c1=1.;c2=0. #Case (a) if l==1: q1=1.;q2=0.;c1=1.;c2=0. #Case (b) if l==2: q1=0.;q2=1.;c1=0.;c2=1. #Case (c) ←(b)は非可制御 (q2=0) ←(c)は非可観測 (c1=0) DGC_no1/ page47v1

3.4 応答の算定（プログラム実行結果） 42 書籍の図 page45 “高橋安人,ディジタル制御 ,岩波書店,1987

3.5 安定に関する定理（リアプノフ式を解く） 3.1 状態空間系 3.2 伝達関数 3.3 同伴形（正準形） 3.4 応答の算定 3.5 安定に関する定理 43

3.5 安定に関する定理（リアプノフ式を解く） 44 リアプノフは状態空間における状態の運動に注目して安定を検討した。彼の定理を線形離散時間系に適用す ると次の式が出てくる。一般にWは単位行列Iまたは正の重みwiをi=1, …, nの対角要素とする対角形をとる。 リアプノフ式： 𝑷𝑇 𝒀𝑷 − 𝒀 = −𝑾 (3-37) リアプノフ式を解くプログラム実行結果 例１ P= [[ 0.6 -0.4 0.3] [ 0.7 -0.9 0.2] [-0.7 0.5 -0.3]] W= [[-1.] [-0.] [-2.] [-0.] [-0.] [-3.]] page53 プログラムでは列ベ クトルで入力する 例2 P= [[ 0.6 -0.4 0.3] [ 0.7 -0.9 0.2] [-0.7 0.5 -0.3]] W= [[-1.] [-0.] [-1.] [-0.] [-0.] [-1.]] 注1) 連続時間系では以下がリアプノフ式 となる 𝑑𝒙(𝑡)Τ𝑑𝑡 = 𝑨𝒙(𝑡) 𝑨𝑇 𝒀 + 𝒀𝑨 = −𝑾 単位行列にしてみた Y= [[ 4.50411 -3.81373 1.14657] [-3.81373 6.8467 -1.02509] [ 1.14657 -1.02509 3.44859]] Y= [[ 2.68074 -1.83847 0.55835] [-1.83847 3.3489 -0.49822] [ 0.55835 -0.49822 1.22405]] 検算 P'YP-Y=-W 検算 P'YP-Y=-W W= [[ 1. -0. -0.] [ 0. 2. -0.] [-0. -0. 3.]] W= [[ 1. -0. 0.] [-0. 1. -0.] [ 0. -0. 1.]] 注2) (3-37)の導出は下記にある。 システムと制御（高橋安人著、 岩波書店、1978，p.207） （後に10章で必要のためここでプログラムを 作成とのこと） DGC_no1/ page50

付録：slycotのInstall 45 これまでにPython-Controlで示した機能を使うには”slycot”Libraryが必要です 詳しくは次のHPで確認ください： https://pypi.org/project/slycot/ 私は”Anaconda”環境でpythonを利用しています。 その場合のinstallは次のようにしました。 1. conda update –all 私の場合は必要でした（不要かもしれませんが） 2. conda install -c conda-forge slycot 3. pytest --pyargs slycot これはキチンと組込めたかのテストです （前提：controlはinstall済） これらのコマンドは、全て Anaconda Prompt画面上で実行しました

• https://pypi.org/project/slycot/

付録：級数展開 sympyを用いて級数展開 1 = 0 + 𝑧 −1 + 𝑝𝑧 −2 + 𝑝 2 𝑧 −3 + 𝑝3 𝑧 −4 , ⋯ 𝑧−𝑝 0, 1, 𝑝 , 𝑝2 , 𝑝3, ⋯ 1 𝑧−𝑝 2 = 0 + 0 ∙ 𝑧 −1 + 𝑧 −2 + 2𝑝𝑧 −3 + 3𝑝2 𝑧 −4 + 4𝑝 3 𝑧 −5 , ⋯ 0, 0, 1, 2𝑝 , 3𝑝2 , 4𝑝3 , ⋯ 1 𝑧 −1 = 𝑧 − 𝑝 1 − 𝑝 ∙ 𝑧 −1 𝑥 = 𝑧 −1 として級数展開した この方法が正しいかは？ 46

つづく ここまで３章です。次回は４章以降に進みます。引き続きプログラ ムのPythonでの再記述を続ける予定です。 • 今回Pythonで用いた変数などは、書籍とできるだけ同じ名前を用いました。このた め、各章プログラム間で変数名の混乱が生じ、プログラムの判読が困難になってき ました。 ✓ これは原著が1980年代のHP-Basicで記述であり、当時は変数名の制限、メモ リ容量の制限などがあり、止むを得ない点であると思います。 ✓ 但し、今回pythonでの再記述では、一部変数名を変更しました。 ✓ pythonは、現在使いながら学んでいます。このため、おかしな点も多々あるか と思います。 • また、計算結果の検証のため、python-controlでの解析結果を示しました。 python-controlは、使い方も簡単で、いろいろなことができることに非常に感心致 しました。 47