ディジタル制御（1987年）その２v2 pythonで記述してみる（４.連続時間プラントの離散時間形）

2024.2.7 ディジタル制御（1987年）その２v2 pythonで記述してみる（４.連続時間プラントの離散時間形） 下 条 誠 電気通信大学名誉教授 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 最適フィルタ

コメント 本書の面白いところは、全てのアルゴリズムが、HP-Basi言語と付属の基本行列 演算で作られている点です。正準形変換の方法やリカチ行列の解法などの記述が 多々あり大変興味深いものです。これらは現在では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/

４. 4.1 連続時間プラントの離散時間形 プラント微分方程式と等価の差分式 4.2 パルス伝達関数と拡張ｚ変換 4.3 解析的手法によるｚ変換と拡張ｚ変換 4.4 時系列からの伝達関数 4.5 観測器による状態推定 4章 連続時間プラント 離散時間形 G(s) G(z) 変換方法の話し 4

4.1 プラント微分方程式と等価の差分式 変 換 必 要 5 𝑮 𝑠 𝑮 𝑧 微分式 𝑑 𝒙 𝑡 = 𝑨𝒙 𝑡 + 𝒃𝑢(𝑡) 𝑑𝑡 問題： 差分式 (4-1) 𝒙 𝑘 + 1 = 𝑷𝒙 𝑘 + 𝒒𝑢(𝑘) A行列(n×n) P行列(n×n) bベクトル(n) qベクトル(n) Tサンプリング周期 変換する 1入出力系 (4-2)

4.1 プラント微分方程式と等価の差分式 6 １）出力は自由応答と強制応答の和となる 入力 𝑢(𝑡) 𝑘T 1： 連続時間 プラント 𝑇 𝑘+1 T 𝒙 𝑘 = 𝑒 𝑘𝑨𝑇 𝒙 0 = 𝑷𝒌 𝒙 0 1. その時の自由応答 2. そして、強制応答 （初期条件成分のなごり） 等比系列成分 2： （時刻ｋ成分のなごり+強制応答） 𝒙 𝑘 + 1 = 𝑷𝒙 𝑘 + 𝒒𝑢(𝑘) 時刻kからの等比成分 強制応答 𝑘+1 𝒙 𝑘+1 = 𝑒 𝑨𝑇 𝑒𝑨 𝒙 𝑘 +න 𝑘𝑇 𝑘+1 𝑇−𝜏 𝒃𝑢 𝜏 𝑑𝜏 (変数置換)注 𝑇 𝒙 𝑘+1 = 𝑒 𝑨𝑇 𝒙 𝑘 + න 𝑒 𝑨𝜂 𝑑𝜂𝒃 𝑢 𝑘 （Aとbは微分式(4-1)の値） 0 P q P, qが分かればよい 注) 𝜂 = 𝑘+1 𝑇−𝜏

4.1 プラント微分方程式と等価の差分式 7 応答を求めるには、Pとq（Q0）が分かればよい 𝑇 𝑷 = 𝑒 𝑨𝑇 𝑸0 = න 𝑒 𝑨𝜂 𝑑𝜂 𝒒 = 𝑸0 𝒃 (4-3) 0 算出方法：行列指数関数の展開形を用いて求める 𝑷 = 𝑒 𝑨𝑇 ≑ 𝐼 + 𝑨𝑇 + 1 𝑨𝑇 2! 𝑇 𝑸0 = න 𝑒 𝑨𝜂 𝑑𝜂 ≑ 𝑻 𝐼 + 0 2 + ⋯+ 1 𝑨𝑇 𝑁! 1 1 𝑨𝑇 + 𝑨𝑇 2! 3! 2 𝑁 + ⋯+ 1 𝑨𝑇 𝑁+1 ! 𝑁 これらの無限数列の計算は、list4-1で(AT)N/N!の大きさが所定の 最小値より小さくなったところで打ち切る DGC_no2/p57v4

4.1 プラント微分方程式と等価の差分式（例題） 微分式 差分式 𝑑 𝒙 𝑡 = 𝑨𝒙 𝑡 + 𝒃𝑢(𝑡) 𝑑𝑡 𝒙 𝑘 + 1 = 𝑷𝒙 𝑘 + 𝒒𝑢(𝑘) (4-1) 次の3連振動系を取り上げる。 例題： この微分式を求め、差分式に変換する y1 (t) u1 (t) y2 (t) b1 m1 8 k1 u2 (t) y3 (t) b2 m2 k2 u3 (t) b3 m3 k3 (4-2)

4.1 プラント微分方程式と等価の差分式（例題） １）3連振動系の微分式を求める y1 (t) u1 (t) m1 page 239 y2 (t) b1 u2 (t) k1 x1=y1=m1の変位 y3 (t) b2 m2 9 u3 (t) k2 x2=dx1/dtその速度 b3 m3 x3=y3=m2の変位 k3 x4=dx3/dtその速度 x5=y3=m1の変位 x6=dx5/dtその速度 𝑑𝑥1 𝑡 Τ𝑑𝑡 = 𝑥2 𝑡 状態式と出力式： 𝑑 𝒙 𝑡 = 𝑨𝒙 𝑡 + 𝑩𝒖(𝑡) 𝑑𝑡 𝒚 𝑡 = 𝑪𝒙 𝑡 𝑚1 𝑑𝑥2 𝑡 Τ𝑑𝑡 = −𝑘1 𝑥1 𝑡 − 𝑏1 𝑥2 𝑡 +𝑘1 𝑥3 𝑡 +𝑏1 𝑥4 𝑡 + 𝑢1 (𝑡) 𝑑𝑥3 𝑡 Τ𝑑𝑡 = 𝑥4 𝑡 𝑚2 𝑑𝑥4 𝑡 Τ𝑑𝑡 = 𝑘1 𝑥1 𝑡 + 𝑏1 𝑥2 𝑡 − 𝑘1 + 𝑘2 𝑥3 𝑡 − 𝑏1 + 𝑏2 𝑥4 𝑡 +𝑘2 𝑥5 𝑡 +𝑏2 𝑥6 𝑡 + 𝑢2 (𝑡) 𝑑𝑥5 𝑡 Τ𝑑𝑡 = 𝑥6 𝑡 𝑚3 𝑑𝑥6 𝑡 Τ𝑑𝑡 = 𝑘3 𝑥3 𝑡 + 𝑏2 𝑥4 𝑡 − 𝑘2 + 𝑘3 𝑥5 𝑡 − 𝑏2 + 𝑏3 𝑥6 𝑡 +𝑢2 (𝑡) DGC_no2/ page239

4.1 プラント微分方程式と等価の差分式（例題） 10 ２）3連振動系の微分式を求め、状態式と出力式を算出する page 239 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 Bとして、外力を質量m1 （自由端）とm3に加える 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 m1=m2=m3=1 とする Cは2出力としてm1 とm2の変位にした （B,Cの数値は書籍のまま） 3連振動系でバネ定数を(k1=1,k2=2,k3=1)、減衰係数を(b1=0.01,b2=0,b3=0.01)とした A matrix [[ 0.000 1.000 0.000 0.000 0.000 0.000] [-1.000 -0.010 1.000 0.010 0.000 0.000] [ 0.000 0.000 0.000 1.000 0.000 0.000] [ 1.000 0.010 -3.000 -0.010 2.000 0.000] [ 0.000 0.000 0.000 0.000 0.000 1.000] [ 0.000 0.000 2.000 0.000 -3.000 -0.010]] B matrix [[0.000 0.000] [1.000 0.000] [0.000 0.000] [0.000 0.000] [0.000 0.000] [0.000 2.000]] C matrix [[1.000 0.000 0.000 0.000 0.000 0.000] [0.000 0.000 1.200 0.000 0.000 0.000]]

4.1 プラント微分方程式と等価の差分式（例題） ２）3連振動系の微分式を求め、状態式と出力式を算出する(python control) 試しに伝達関数Gij(s)を表示させてみる 伝達関数 11

4.1 プラント微分方程式と等価の差分式（例題） 12 ３）3連振動系の微分式を求め、状態式と出力式を算出する（List4-1） 微分方程式と等価の差分式への変換プログラム実行結果 -------------A, B and C matrix------------ A matrix [[ 0.000 1.000 0.000 0.000 0.000 0.000] [-1.000 -0.010 1.000 0.010 0.000 0.000] [ 0.000 0.000 0.000 1.000 0.000 0.000] [ 1.000 0.010 -3.000 -0.010 2.000 0.000] [ 0.000 0.000 0.000 0.000 0.000 1.000] [ 0.000 0.000 2.000 0.000 -3.000 -0.010]] B matrix [[0.000 0.000] [1.000 0.000] [0.000 0.000] [0.000 0.000] [0.000 0.000] [0.000 2.000]] C matrix [[1.000 0.000 0.000 0.000 0.000 0.000] [0.000 0.000 1.200 0.000 0.000 0.000]] 𝑑 𝒙 𝑡 = 𝑨𝒙 𝑡 + 𝑩𝑢(𝑡) 𝑦 𝑡 = 𝑪𝒙 𝑡 𝑑𝑡 -------------P, Q and dc-gain matrix-----------number of recursion= 11 P matrix [[ 0.88047 [-0.45762 [ 0.11458 [ 0.41907 [ 0.00490 [ 0.03800 0.47853 0.11420 0.87589 0.41688 0.02092 0.66071 0.11878 -1.22197 0.00055 0.21976 0.00528 0.76430 0.02092 0.00528 0.11878 0.04019 0.43999 0.21976 0.65652 0.76435 0.03854 0.65808 0.21938 -1.24198 0.00055] 0.00528] 0.03854] 0.21938] 0.43969] 0.65368]] Q matrix [[0.12224 0.00009] [0.47853 0.00109] [0.00271 0.00989] [0.02092 0.07709] [0.00005 0.23454] [0.00055 0.87938]] dc-gain [[2.50000 2.00000] [1.80000 2.40000]] 𝒙 𝑘 + 1 = 𝑷𝒙 𝑘 + 𝒒𝑢(𝑘) （List4-1） DGC_no2/p57v4

4.1 プラント微分方程式と等価の差分式（例題） ４）3連振動系の微分式を求め、状態式と出力式を算出する（PLOT） 13 page 59 （2入力として力をm1とm3に加え、2出力としてm1とm2の変位とした） m3への入力 による変位 m1への入力 による変位 m3への入力 による変位 m1への入力 による変位 DGC_no2/ p59v2

4.1 プラント微分方程式と等価の差分式（例題） 14 書籍の図 page57 “高橋安人,ディジタル制御 ,岩波書店,1987

4.1 プラント微分方程式と等価の差分式(確認) 15 ５）3連振動系の微分式を求め、状態式と出力式を算出する（確認:python-control） import control as ct Tsample=0.5 sys=ct.ss(A,B,C,0) #print(sys) sysPd=ct.c2d(sys, Tsample, method='zoh') print(sysPd) 連続形から離散形への変換 -------------P, Q and dc-gain matrix-----------number of recursion= 11 A = [[ 0.88047 [-0.45762 [ 0.11458 [ 0.41907 [ 0.00490 [ 0.03800 P B = [[0.12224 [0.47853 [0.00271 [0.02092 [0.00005 [0.00055 Q 0.47853 0.11420 0.87589 0.41688 0.02092 0.66071 0.11878 -1.22197 0.00055 0.21976 0.00528 0.76430 0.00009] 0.00109] 0.00989] 0.07709] 0.23454] 0.87938]] 0.02092 0.00528 0.11878 0.04019 0.43999 0.21976 0.65652 0.76435 0.03854 0.65808 0.21938 -1.24198 0.00055] 0.00528] 0.03854] 0.21938] 0.43969] 0.65368]] 上記(control)を用いて 計算したもの C = [[1.00000 0.00000 0.00000 0.00000 0.00000 0.00000] [0.00000 0.00000 1.20000 0.00000 0.00000 0.00000]] D = [[0.00000 0.00000] [0.00000 0.00000]] P matrix [[ 0.88047 [-0.45762 [ 0.11458 [ 0.41907 [ 0.00490 [ 0.03800 0.47853 0.11420 0.87589 0.41688 0.02092 0.66071 0.11878 -1.22197 0.00055 0.21976 0.00528 0.76430 Q matrix [[0.12224 0.00009] [0.47853 0.00109] [0.00271 0.00989] [0.02092 0.07709] [0.00005 0.23454] [0.00055 0.87938]] dc-gain [[2.50000 2.00000] [1.80000 2.40000]] 0.02092 0.00528 0.11878 0.04019 0.43999 0.21976 0.65652 0.76435 0.03854 0.65808 0.21938 -1.24198 0.00055] 0.00528] 0.03854] 0.21938] 0.43969] 0.65368]] 本章の方法で計算 したもの DGC_no2/p57v4

4.1 プラント微分方程式と等価の差分式(確認) 実行：Jupyter Notebook import control as ct Tsample=0.5 連続形 16 離散形状態空間モデルに変更する sys=ct.ss(A,B,C,0) #print(sys) 離散形 sysPd=ct.c2d(sys, Tsample, method='zoh');sysPd MIMO 記述子状態空間モデル A B C D 試しに伝達関数Gij(z)を表示させてみる G(z) 多入出力系の伝達関数Gij分母は全て同じ

• https://jp.mathworks.com/help/control/ug/mimo-state-space-models.html#bsst_fp-1

4.2 パルス伝達関数と拡張ｚ変換 4.1 プラント微分方程式と等価の差分式 4.2 パルス伝達関数と拡張ｚ変換 4.3 解析的手法によるｚ変換と拡張ｚ変換 4.4 時系列からの伝達関数 4.5 観測器による状態推定 𝑮 𝑠 𝑮 𝑧 数値計算で求める 17

4.2 パルス伝達関数と拡張ｚ変換 18 ホールドとサンプラのついた連続時間プラントのG(s)から、前スライド図4-1(b)のパルス 伝達関数G(z)を求める方法は、以下の2つがある。ここでは１の方法を述べる 1. 状態空間における数値計算 2. 解析的手法 𝑮 𝑠 = 𝒀(𝑠) 𝑼(𝑠) 𝑮 𝑧 = 求め方 Step1: 𝑑 𝒙 𝑡 = 𝑨𝒙 𝑡 + 𝒃𝑢(𝑡) 𝑑𝑡 𝒚 𝑡 = 𝒄𝒙 𝑡 Step2: 𝑨, 𝒃 (list 4-1) 注 𝑷, 𝒒 ( 4-1) 微分式から差分式へ変換 p57v1.py Step3: 𝑷, 𝒒, 𝒄 𝒀(𝑧) 𝑼(𝑧) 𝑷𝒄 , 𝒒𝒄 , 𝒄𝒄 (list 3-1) page40v1.py 同伴形を導出する 正準形のこと 注）多入出力系の同伴形はそう簡単ではない。本章では1入出力系とする

4.2 パルス伝達関数と拡張ｚ変換 「遅れ時間なし」と「遅れ時間あり」に分けて考える case (a) L1=0 case (b) 0<L1<T y(k+1) y(k) u(k) u(k) u(k-1) u(k-1) kT y(k+1) y(k) L1 (k+1)T t (k-1)T (k-1)T (k+1)T t kT 𝑄2 𝑸0 算出区間 T-L1 𝑄1 L1 :時間遅れ 連続時間プラントの入出力間にむだ 時間遅れL1があると、階段状入力が 上図のようにズレる。 伝達関数はQ0から決まる 伝達関数はQ1とQ2の和から決まる 19

4.2 パルス伝達関数と拡張ｚ変換（遅れ時間なし） case (a) (a)L1=0 L1=0 20 (b) 0<L1<T y(k+1) y(k) y(k+1) y(k) u(k) u(k) u(k-1) u(k-1) kT L1 (k+1)T t (k-1)T (k-1)T kT T-L1 (k+1)T t 𝑸0 算出区間 𝑑 𝒙 𝑡 = 𝑨𝒙 𝑡 + 𝒃𝑢(𝑡) 𝑑𝑡 微分式 従 来 通 り (list 4-1) 差分式 𝑇 𝒒 = න 𝑒 𝑨𝜂 𝑑𝜂 𝒃 = 𝑸0 𝒃 page57v4.py 0 𝒙 𝑘 + 1 = 𝑷𝒙 𝑘 + 𝒒𝑢(𝑘) 𝑷𝒄 = 𝑻𝑷𝑻−1 𝒒𝒄 = 𝑻𝒒 𝑪𝒄 = 𝑪𝑻−1 𝑪𝒄 = 𝑏𝑛 0 0 0 𝒒𝒄 = ⋮ 0 1 𝑏𝑛−1 𝑏𝑛−2 (list 3-1) page40v1.py ｚ変換式 同伴形（正準形）に変換 0 0 0 𝑷𝒄 = ⋮ 0 −𝑎𝑛 1 0 0 ⋮ 0 −𝑎𝑛−1 0 1 0 ⋮ 0 ⋯ … 0 ⋯ 0 ⋯ 0 ⋮ ⋮ ⋯ 1 ⋯ −𝑎1 ⋯ 𝑏1

4.2 パルス伝達関数と拡張ｚ変換（遅れ時間あり） 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 𝑒 𝑨𝜂 𝑑𝜂 21 𝑇 𝑒 𝑨𝜂 𝑑𝜂 = 𝑸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) 𝑎

4.2 パルス伝達関数と拡張ｚ変換（遅れ時間あり) 微分式： 差分式： 𝑑 𝒙 𝑡 = 𝑨𝒙 𝑡 + 𝒃𝑢(𝑡) 𝑑𝑡 22 y(k+1) y(k) u(k) 𝒙 𝑘 + 1 = 𝑷𝒙 𝑘 + 𝒒𝑢(𝑘) 自由応答項 𝑒 𝐴𝑇 𝒙 𝑘 u(k-1) L1 強制応答項 𝑸0 𝒃 (k-1)T u(k-1)で の応答 0 (k+1)T t kT 𝑇 න 𝑒 𝐴𝜂 𝑑𝜂 ∙ 𝒃 𝑢(𝑘) T-L1 𝑄2 𝑄1 u(k)での 応答 𝑸0 = 𝑸1 + 𝑸2 𝒒 以上から次のようにも表せる 差分式： 𝒙 𝑘 + 1 = 𝑷𝒙 𝑘 + 𝒒1 𝑢 𝑘 + 𝒒2 𝑢(𝑘 − 1) 1個時間遅れ入力 伝達関数： 伝達関数はQ1と Q2の和になる 𝐵1 𝑧 𝐵2 𝑧 −1 𝑧𝐵1 𝑧 + 𝐵2 𝑧 𝐺 𝑧 = + 𝑧 = 𝐴 𝑧 𝐴 𝑧 𝑧𝐴 𝑧 𝑄1 𝑄2 1個時間遅れ入力 次数が n-->n+1 となる

4.2 パルス伝達関数と拡張ｚ変換（例題) 23 対象システム：3連振動系（n=6のシステム） y1 (t) u1 (t) y2 (t) u2 (t) b1 y3 (t) b2 u3 (t) b3 算出例を次頁以降に示す m1 m2 k1 k2 m3 k3 3連振動系の状態式 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 𝑘1 0 − 𝑘1 + 𝑘2 0 𝑘2 0 𝑏1 1 − 𝑏1 + 𝑏2 0 𝑏2 0 0 0 𝑘2 0 − 𝑘2 + 𝑘3 0 0 0 𝑏2 1 − 𝑏2 + 𝑏3 0 1.0 0 B= 1.0 0 1.0 C = 1. 1. 1. 1. 1. 1.

4.2 パルス伝達関数と拡張ｚ変換（遅れ時間なし） DGC_no2/ p63q1q2 output 確認OK P= [[ 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]] q1= q2= [[0.125 ] [[0.] [0.5 ] [0.] [0.125 ] [0.] [0.4994] [0.] [0.1223] [0.] [0.4788]] [0.]] A1(z):(an,an-1...a2,a1:) [ 0.9851 -4.3394 9.0348 -11.3515 B1(z):(bn..b1:) [-1.0915 4.536 -8.2534 8.95 9.0831 -4.3854] -5.6703 1.8504] ----- G(z)=(zB1(z)+B2(z))/zA(z) ----B(z)= [ 0. -1.0915 4.536 -8.2534 8.95 -5.6703 1.8504] 書籍表4-3とも同値 polynominal:(分子） 6 5 4 3 2 1.85 z - 5.67 z + 8.95 z - 8.253 z + 4.536 z - 1.092 z /zA(z) python-control output 24 Tsample=0.5

4.2 パルス伝達関数と拡張ｚ変換（時間遅れあり） 今回の計算結果は書籍の値と違う値となりました。 何回か見直しましたが、その原因が分からないため、 計算結果をそのまま以下に示します。 但し、遅延時間ゼロ（L1=0)での値は一致しています 25

書籍の表 26 書籍の表 ここは同じ これ以外は 値が違う “高橋安人,ディジタル制御 ,岩波書店,1987 page64

4.2 パルス伝達関数と拡張ｚ変換：計算結果 L1=0.0 q1= [[0.125 ] [0.5 ] [0.125 ] [0.4994] [0.1223] [0.4788]] L1=0.2 q2= [[0.] [0.] [0.] [0.] [0.] [0.]] zB1(z)+B2(z): 6 5 4 3 2 1.85 z - 5.67 z + 8.95 z - 8.253 z + 4.536 z - 1.092 z q1= [[0.045 ] [0.3 ] [0.045 ] [0.3 ] [0.0446] [0.2951]] 27 G(z)=(zB1(z)+B2(z))/zA(z) q2= [[0.08 ] [0.2 ] [0.08 ] [0.1995] [0.0776] [0.1837]] zB1(z)+B2(z): 6 5 4 3 2 1.03 z - 2.288 z + 2.303 z - 0.3473 z - 1.413 z + 1.568 z - 0.5309 （分子、分母からｚが相殺されることになる） L1=0.1 q1= [[0.08 ] [0.4 ] [0.08 ] [0.3998] [0.0788] [0.3888]] zB1(z)+B2(z): L1=0.3 q2= [[0.045 ] [0.1 ] [0.045 ] [0.0996] [0.0434] [0.09 ]] 6 5 4 3 2 1.427 z - 3.921 z + 5.5 z - 4.138 z + 1.429 z + 0.3039 z - 0.2806 Tsample=0.5 q1= [[0.02 ] [0.2 ] [0.02 ] [0.2 ] [0.0199] [0.1985]] q2= [[0.105 ] [0.3 ] [0.105 ] [0.2995] [0.1023] [0.2803]] zB1(z)+B2(z): 6 5 4 3 2 0.6584 z - 0.7791 z - 0.6313 z + 3.105 z - 3.98 z + 2.698 z - 0.7499 DGC_no2/ p63q1q2

4.2 パルス伝達関数と拡張ｚ変換：計算結果 ゼロ点の計算結果：書籍とprogramの値の比較 DGC_no2/ p63q1q2 28 Tsample=0.5 書籍のようにL1=0.1~0.3で負の 実数ゼロが大きくは変化しない DGC_no2/ page63_poly

4.2 パルス伝達関数と拡張ｚ変換（遅れ時間なし） 計算した伝達関数を用いて出力結果をPLOTした 伝達関数： 𝐺 𝑧 = Tsample=0.5 𝐵1 𝑧 𝐵2 𝑧 −1 + 𝑧 𝐴 𝑧 𝐴 𝑧 𝑄1 遅れ時間なし L1=0 29 𝑄2 1個時間遅れ入力 impulse 応答 impulse input DGC_no2/ p63resp

4.2 パルス伝達関数と拡張ｚ変換（遅れ時間あり） 30 Tsample=0.5 遅れ時間あり L1≠0 impulse L1=0.1 L1=0.3 L1=0.2 L1=0.4 impulse 応答 DGC_no2/ p63resp

4.3 4.1 解析的手法によるｚ変換と拡張ｚ変換 プラント微分方程式と等価の差分式 4.2 パルス伝達関数と拡張ｚ変換 4.3 解析的手法によるｚ変換と拡張ｚ変換 4.4 時系列からの伝達関数 4.5 観測器による状態推定 𝑮 𝑠 𝑮 𝑧 1つの方法について記述 1. 解析的手法 2. 実測値 (g1,g2,…gn) からのG(z)推定 3. 実測値 (g1,g2,…gn) からのP, q推定 31

4.3 解析的手法によるｚ変換と拡張ｚ変換 1. 状態空間における数値計算 2. 解析的手法 ←４.２ ←４.３ 微分式 差分式 𝑑 𝒙 𝑡 = 𝑨𝒙 𝑡 + 𝑩𝑢(𝑡) 𝑑𝑡 𝒚 𝑡 = 𝑪𝒙 𝑡 + 𝑫𝑢(𝑡) list 4-1 𝒙 𝑘 + 1 = 𝑷𝒙 𝑘 + 𝑸𝑢(𝑘) 𝒚 𝑘 = 𝑪𝒙 𝑘 + 𝑫𝑢(𝑘) 正準形のこと page57v4.py 𝓏 ℒ ｚ変換式 Laplace変換式 𝐺 𝑧 = 1−𝑧 𝒀(𝑠) 𝑮 𝑠 = 𝑼(𝑠) 32 −1 𝒵 ℒ −1 1 𝐺 𝑠 𝑠 同伴形 list 3-1 page40v1.py 多項式伝達関数 １入出力型 (4-8) 𝒀(𝑧) 𝑮 𝑧 = 𝑼(𝑧) （連続型のs関数を、離散型のz関数に近似する方法として双一次変換もある）

4.3 解析的手法によるｚ変換と拡張ｚ変換 G(s)とG(z)の間には右の関係が成立する。①1/s G(s)は 𝐺 𝑧 = 1 − 𝑧 −1 𝒵 ℒ −1 連続時間形での単位ステップ応答であり。それを②離散 33 1 𝐺 𝑆 𝑆 ① 時間形に変換して、③一刻みの時間差をとることで、単 連続時間形ステップ応答 ② 位インパルス応答としている 離散時間形ステップ応答 ③ 離散時間形インパルス応答 page66 “高橋安人,ディジタル制御 ,岩波書店,1987

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

4.3 解析的手法によるｚ変換と拡張ｚ変換（例） 1. 解析的手法 𝐾𝑒 −𝑠𝐿 𝐺 𝑠 = 1 + 𝑇1 𝑠 1 + 𝑇2 𝑠 どちらも遅れ時間L=3の応答 (T:sampling interval) ただし、サンプリング時間がT=2, T=4の違いがある T=2, N=1, L1=1 (サンプリング、1個遅れ+L1) T=4, N=0, L1=3 横軸はstep数なので時間軸にすると2倍違う 注意）g(k)は5倍して表示した DGC_no2/ p70T2_g(z) DGC_no2/ p70T4_g(z) 35

4.4 時系列からの伝達関数 4.1 プラント微分方程式と等価の差分式 4.2 パルス伝達関数と拡張ｚ変換 4.3 解析的手法によるｚ変換と拡張ｚ変換 4.4 時系列からの伝達関数 4.5 観測器による状態推定 𝑮 𝑠 𝑮 𝑧 2つの方法について記述 1. 解析的手法 2. 実測値 (g1,g2,…gn) からのG(z)推定 3. 実測値 (g1,g2,…gn) からのP, q推定 36

4.4 時系列からの伝達関数 37 連続時間プラントの単位ステップ応答の実測値からのG(z)推定する。 ∞ 𝐺 𝑧 = 𝑔1 𝑧 −1 + 𝑔2 𝑧 −2 + ⋯ = ෍ 𝑔𝑗 𝑧 −𝑗 𝑔𝑗 = 𝑦𝑗 − 𝑦𝑗−1 𝑗 = 1,2, ⋯ (4-13) 𝐽=1 nとpを決める 𝐺 𝑧 −1 𝑝1 = 𝑒= 𝑏1 𝑧 −1 + 𝑏2 𝑧 −2 + ⋯ + 𝑏𝑛 = 1 − 𝑝𝑧 −1 𝐺 1 − 𝑔1 + 𝑔2 + ⋯ + 𝑔𝑛 𝐺 1 − 𝑔1 + 𝑔2 + ⋯ + 𝑔𝑛−1 𝑝1 − 𝑝2 𝑝2 < 0.05 𝑏1 = 𝑔1 , 𝑏2 = 𝑔2 − 𝑝𝑔1 , ⋯ , 𝑏𝑛 = 𝑔𝑛 − 𝑝𝑔𝑛−1 𝑝2 = 𝑔𝑛+1 𝑔𝑛 程度となるnとpとする 次頁に4.3で用いた例で計算したものを示す nとpを決める方法 (4-19)

4.4 時系列からの伝達関数 ２．連続時間プラントの単位ステップ応答の実測値 (g1,g2,…gn) からのG(z)推定 𝐺 𝑧 −1 𝑏1 𝑧 −1 + 𝑏2 𝑧 −2 + ⋯ + 𝑏𝑛 = 1 − 𝑝𝑧 −1 𝑏1 = 𝑔1 , 𝑏2 = 𝑔2 − 𝑝𝑔1 , ⋯ , 𝑏𝑛 = 𝑔𝑛 − 𝑝𝑔𝑛−1 横軸はstep数なので時間軸にすると2倍違う DGC_no2/ p71T2_gi DGC_no2/ p71T4_gi 38

時系列からの状態方程式形 39 時系列から伝達関数、そして状態空間形への変換 例）状態空間形への変換（1入力系） 伝達関数： 𝑏1 𝑧 𝑛−1 + 𝑏2 𝑧 𝑛−2 + ⋯ + 𝑏𝑛 G 𝑧 = 𝑧 𝑛−1 𝑧 − 𝑝 𝑏1 𝑧 −1 + 𝑏2 𝑧 −2 + ⋯ + 𝑏𝑛 −1 𝐺 𝑧 = 1 − 𝑝𝑧 −1 (4-19) 𝑏1 = 𝑔1 , 𝑏2 = 𝑔2 − 𝑝𝑔1 , ⋯ , 𝑏 𝑛 = 𝑔𝑛 − 𝑝𝑔𝑛−1 状態方程式： 0 1 0 0 ⋮ ⋯ 𝑃= 0 0 0 0 0 0 𝒙 𝑘 + 1 = 𝑷𝒙 𝑘 + 𝒒𝑢(𝑘) 0 1 ⋯ 0 0 0 0 0 0 0 ⋯ ⋯ 0 1 0 0 0 0 0 𝑔1 0 𝑔2 ⋮ 𝑔3 𝑞 = 0 ⋮ 1 ⋮ 𝑝 𝑔𝑛 𝑦 𝑘 = 𝑪𝒙 𝑘 𝐶= 1 0 0 ⋯ ⋯ 0

時系列からの状態方程式形（例） 40 ３．連続時間プラントの単位ステップ応答の実測値 (g1,g2,…gn) からのP, q推定 0 1 0 0 ⋮ ⋯ 𝑃= 0 0 0 0 0 0 0 0 0 0 1 0 0 0 ⋯ ⋯ ⋯ ⋮ 0 0 1 0 0 0 0 1 0 0 0 𝑝 𝑔1 𝑔2 𝑔3 𝑞= ⋮ ⋮ 𝑔𝑛 𝐶= 1 0 0 ⋯ ⋯ 0 横軸はstep数なので時間軸にすると2倍違う DGC_no2/ p72T2_Pq DGC_no2/ p72T4_Pq

時系列からの状態方程式形（多入出力系） 前頁は、n>mの多入出力系へ拡張可能 𝑩1 𝑧 𝑛−1 + 𝑩2 𝑧 𝑛−2 + ⋯ + 𝑩𝑛 𝑮 𝑧 = ∈ ℝ𝑚×𝑚 𝑛−1 𝑧 𝑧𝑰 − 𝑹 𝟎 𝟎 ⋮ 𝑷= 𝟎 𝟎 𝟎 𝑰 𝟎 𝟎 𝑰 ⋯ ⋯ 𝟎 𝟎 𝟎 𝟎 𝟎 𝟎 𝟎 𝟎 ⋯ 𝟎 𝟎 𝟎 𝟎 𝟎 ⋯ 𝑰 𝟎 𝟎 𝟎 𝟎 ⋮ 𝟎 𝑮𝑛 𝑹 𝑮1 𝑮2 𝑮 𝑸= 3 ⋮ ⋮ 𝑮𝑛 𝑩1 = 𝑮1 , 𝑩2 = 𝑮2 − 𝑮1 𝑹, ⋯ , 𝑩𝑛 = 𝑮𝑛 − 𝑮𝑛−1 𝑹 𝑪= 𝑰 𝟎 𝟎 ⋯ ⋯ 𝟎 (4-22) 𝑷 ∈ ℝ𝑚𝑛×𝑚𝑛 , 𝑸 ∈ ℝ𝑚𝑛×𝑚 , 𝑪 ∈ ℝ𝑚×𝑚𝑛 例）m=2 𝑹= 𝑝1 0 0 𝑝2 𝑮𝒊 = 𝑔11 (𝑖) 𝑔12 (𝑖) 𝑔21 (𝑖) 𝑔22 (𝑖) 𝑖 = 1,2 ⋯ , 𝑛 Rが対角であるのは、応答に含まれる等比成分だけから決めることを意味している。 例えば、u1→y1, u1→y2の応答成分がともに等比p1をもって減衰するものと考える 41

付録：4.3のプラントモデルについて 加熱室とその室内に加熱物体を置き、物体を加熱する例を示す。 室内の温度をx1 (t) 、物体温度をx2 (t)とする 1. 熱収支は 𝑢 𝑡 𝑑𝑥1 𝑡 −𝑥1 𝑡 + 𝑣(𝑡) 𝑥2 𝑡 − 𝑥1 𝑡 = + 𝑑𝑡 𝑅1 𝑅2 𝑑𝑥2 𝑡 𝑥1 𝑡 − 𝑥2 𝑡 物体 ： 𝑐2 = 𝑑𝑡 𝑅2 加熱室： 𝑐1 + 𝑢(𝑡) 𝑋2 (𝑠) 𝑅2 = 𝑋1 (𝑠) 𝑇2 𝑠 + 1 𝐶1 𝑅2 𝐶2 𝑥1 𝑡 𝑥2 𝑡 CR型プラント 2. C1>>C2ならば、加熱室の熱収支で物体の存在を無視できる 𝑋1 (𝑠) 𝑅1 = , 𝑈(𝑠) 𝑇1 𝑠 + 1 𝑅1 𝑣 𝑡 ～ ～ ～ 42 𝑇1 = 𝐶1 𝑅1 𝑇2 = 𝐶2 𝑅2 C1：加熱室の熱容量 R1：加熱室壁の熱伝達抵抗 C2：加熱物体の熱容量 R2：物体表面の熱伝達抵抗 v(t)：外気温度（外乱） page235 3. 加熱物体の温度を出力y(t)とすると、全体の伝達関数は下記となる 𝑌 𝑠 = 𝑋2 (𝑠) 𝐺 𝑠 = 𝑌(𝑠) 𝑅1 𝑅2 ≃ 𝑈(𝑠) 𝑇1 𝑠 + 1 𝑇2 𝑠 + 1 この様な関係は、圧力、拡散、混合、 反応過程にも適用されるそうです

付録：4.3のプラントモデルについて 化学工業プラントの特性はしばしば 次の形の伝達関数で近似される 43 1 + 𝑇1′ 𝑠 1 + 𝑇2′ 𝑠 ⋯ 1 + 𝑇𝑚′ 𝑠 𝐺 𝑠 =K 1 + 𝑇1 𝑠 1 + 𝑇2 𝑠 ⋯ 1 + 𝑇𝑛 𝑠 (n>m) G(s)の単位ステップ応答は0から立上り、 ゲインKへ暫定して頭を打つ（０型） 左図に単位ステップ応答の例を示す。こ の4次系で短い時定数2と1の和3をむだ時 間に置き換え、長い時定数6と10はその ままにすると、次の単位ステップ応答は 4次系のものとほとんど重なる 𝐾𝑒 −3𝑠 𝐺 𝑠 = 1 + 10𝑠 1 + 6𝑠 工業用プロセスではこの原理を使って、 ステップ応答の実験結果に合うように K,L,T1,T2を決めることが良くあるとの こと page236 𝐾𝑒 −𝑠𝐿 𝐺 𝑠 = 1 + 𝑇1 𝑠 1 + 𝑇2 𝑠

4.5 観測器による状態推定 4.1 プラント微分方程式と等価の差分式 4.2 パルス伝達関数と拡張ｚ変換 4.3 解析的手法によるｚ変換と拡張ｚ変換 4.4 時系列からの伝達関数 4.5 観測器による状態推定 44

4.5 観測器による状態推定 プラント： モデル： 𝒙 𝑘 + 1 = 𝑷𝒙 𝑘 + 𝒒𝑢(𝑘) ෝ 𝑘 + 1 = 𝑷ෝ 𝒙 𝒙 𝑘 + 𝒒𝑢(𝑘) 45 一般に 𝒙 0 ≠ 𝟎 ෝ 0 =𝟎 𝒙 推定値誤差の収束を早めるアルゴリズム 暫定推定値 step1: ෝ° 𝑘 + 1 = 𝑷ෝ 𝒙 𝒙 𝑘 + 𝒒𝑢(𝑘) step2: ෝ 𝑘+1 =𝒙 ෝ° 𝑘 + 1 + 𝒇 𝑦 𝑘 + 1 − 𝑦ො ° 𝑘 + 1 𝒙 𝑦ො ° 𝑘 + 1 = 𝒄ෝ 𝒙° 𝑘 + 1 観測器設計問題は重みベクトル f の選定に帰着する。一般にｆは観測誤差ベクトルe(k)の 挙動を支配する。この挙動を観測器系の固有値で指定すると次のようになる 𝒄𝑷 𝒄𝑷2 𝒇 = 𝐹(𝑷) 𝒄𝑷3 ⋮ 𝒄𝑷𝑛 −1 0 0 ⋮ 0 1 (4-25) 𝒄 ∈ ℝ1×𝑛

4.5 観測器による状態推定 46 ここで行列多項式F(P)が指定固有値に関係する。 𝐹 𝑷 = 𝑷𝑛 + 𝛼1 𝑷𝑛−1 + ⋯ + 𝛼𝑛−1 𝑷 + 𝛼𝑛 𝑷 (4-27) において、係数α1…αnを指定固有値を根とする下記の特性式のものとする 𝐹 𝑧 = 𝑧 𝑛 + 𝛼1 𝑧 𝑛−1 + ⋯ + 𝛼𝑛−1 𝑧 + 𝛼𝑛 全ての固有値が0の有限整定観測器ならばF(z)=znだから、下記となる 𝐹(𝑷)＝ 𝑷𝑛 𝒇 = 𝑷𝑛 𝒄𝑷 𝒄𝑷2 𝒄𝑷3 ⋮ 𝒄𝑷𝑛 −1 0 0 ⋮ 0 1 4.4節時系列からの状態方程式形 では行列Pのランクが下がる 逆行列が存在しない？ (4-25)

4.5 観測器による状態推定（例） 47 対象システム：3連振動系（n=6のシステム） y1 (t) u1 (t) 観測器による状態推定を 次頁以降に示す y2 (t) u2 (t) b1 m1 b2 m2 k1 y3 (t) u3 (t) k2 b3 m3 k3 3連振動系の状態式 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 𝑘1 0 − 𝑘1 + 𝑘2 0 𝑘2 0 𝑏1 1 − 𝑏1 + 𝑏2 0 𝑏2 0 0 0 𝑘2 0 − 𝑘2 + 𝑘3 0 0 0 𝑏2 1 − 𝑏2 + 𝑏3 0 1.0 0 B= 0 0 0 C = 1. 0 ステップ入力とした 0 0 0 0

4.5 観測器による状態推定（例） 3連振動系：n=6のシステム 48 # 外乱？としての初期値 X0[0]=np.array([[5],[2],[7],[5],[3],[8]])#書籍表4-6 n=6のシステム 有限整定観測では6刻 みで収束する f : feedback （書籍と一致） DGC_no2/ p77v3

X0[0]=np.array([[5],[2],[7],[5],[3],[8]] 4.5 観測器による状態推定（例） 3連振動系：n=6のシステム 状態値(青)とその 推定値(赤)を重ね て表示 n=6のシステム 有限整定観測では 6刻みで収束する 但し入力が大きい？ f : feedback DGC_no2/ p77v3 49

X0[0]=np.array([[5],[2],[7],[5],[3],[8]] 3連振動系：n=6のシステム feedback=0 としてみた 真値に予測値が 収束しない f : feedback グラフの縦軸が違う ものがある 50

つづく ここまで４章です。次回は5章以降に進みます。 • 今回Pythonで用いた変数などは、書籍とできるだけ同じ名前を用たため、各 章プログラム間で変数名の混乱が一部生じました。 ✓ 今回pythonでの再記述で、変数名を少々変更しました。 • また、書籍にプログラムリストのないものは、記述リストを組合せて私が作 成しました。version1では多くの誤解やミスなどがあり今回訂正しました。 しかし、まだ多くありそうです。 • 前Ver.と比べ、内容は少し詳しく記述しました。 （Pythonのprogramは下記のGitHubに格納してあります） https://github.com/m4881shimojo/YTakahashi_Digital-control/tree/main 51

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

付 録 52

Tsample=0.5 連 続 形 4.2 パルス伝達関数と拡張 ｚ変換（遅れ時間なし） reachable Canonical Form 付録で解説 離 散 形 DGC_no2/p57v4 のシステム 53

4.2 パルス伝達関数と拡張 ｚ変換（遅れ時間なし） Tsample=0.5 連 続 形 離 散 DGC_no2/p57v4 のシステム 形 54

付録：プラント微分方程式と等価の差分式（例題） 55 例) 3連振動系：3入力3出力の例 0 −𝑘1 0 𝐴= 𝑘1 0 0 1 −𝑏1 0 𝑏1 0 0 0 𝑘1 0 − 𝑘1 + 𝑘2 0 𝑘2 0 𝑏1 1 − 𝑏1 + 𝑏2 0 𝑏2 0 0 0 𝑘2 0 − 𝑘2 + 𝑘3 微分方程式と等価の差分式への変換プログラム実行結果 0 0 0 𝑏2 1 − 𝑏2 + 𝑏3 T=0.5 0 1 0 𝐵= 0 0 0 0 0 0 1 0 0 -------------P, Q and dc-gain P matrix [[ 0.88047 0.47853 0.11420 [-0.45762 0.87589 0.41688 [ 0.11458 0.02092 0.66071 [ 0.41907 0.11878 -1.22197 [ 0.00490 0.00055 0.21976 [ 0.03800 0.00528 0.76430 B matrix [[0.000 0.000 0.000] [1.000 0.000 0.000] [0.000 0.000 0.000] [0.000 1.000 0.000] [0.000 0.000 0.000] [0.000 0.000 1.000]] Q matrix [[0.12224 0.00271 0.00005] [0.47853 0.02092 0.00055] [0.00271 0.11730 0.00494] [0.02092 0.43999 0.03854] [0.00005 0.00494 0.11727] [0.00055 0.03854 0.43969]] C matrix [[1.000 0.000 0.000 0.000 0.000 0.000] [0.000 0.000 1.000 0.000 0.000 0.000] [0.000 0.000 0.000 0.000 1.000 0.000]] m1=m2=m3=1 1 𝐶= 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 page 239 k1=1.;k2=2.;k3=1.;b1=0.01;b2=0;b3=0.01 -------------A, B and C matrix-----------A matrix [[ 0.000 1.000 0.000 0.000 0.000 0.000] [-1.000 -0.010 1.000 0.010 0.000 0.000] [ 0.000 0.000 0.000 1.000 0.000 0.000] [ 1.000 0.010 -3.000 -0.010 2.000 0.000] [ 0.000 0.000 0.000 0.000 0.000 1.000] [ 0.000 0.000 2.000 0.000 -3.000 -0.010]] 𝑑 𝒙 𝑡 = 𝑨𝒙 𝑡 + 𝑩𝑢(𝑡) 𝑑𝑡 0 0 0 0 0 1 dc-gain [[2.50000 1.50000 1.00000] [1.50000 1.50000 1.00000] [1.00000 1.00000 1.00000]] matrix-----------0.02092 0.00528 0.11878 0.04019 0.43999 0.21976 0.65652 0.76435 0.03854 0.65808 0.21938 -1.24198 0.00055] 0.00528] 0.03854] 0.21938] 0.43969] 0.65368]] 𝒙 𝑘 + 1 = 𝑷𝒙 𝑘 + 𝒒𝑢(𝑘) DGC_no2/ page239

付録：プラント微分方程式と等価の差分式（例題） 3連振動系： 3入力3出力の例 DGC_no2/ page239 56 外力をm1､m2とm3とに加え、出力をm1､m2とm3 の変位とした このようなシステムでは伝達関数形よりも状態空間形の方が便利である

付録： Reachable Canonical Form？ 57 ChatGPT 「Reachable Canonical Form」（到達可能正準形）は、制御工学や状態空間法において使用される概念の一つ です。状態空間法は、制御システムの挙動をモデル化する手法の一つであり、システムの状態や入力に基づい てシステムの挙動を記述します。 到達可能正準形は、与えられた状態空間表現を特定の形式に変換する手法であり、変換後の形式は特定の性 質を持ちます。この形式に変換されることで、システムの解析や設計が容易になる場合があります。 到達可能正準形への変換の手順にはいくつかのステップがありますが、主な目的は以下の通りです： 1. 到達可能性の確認: • まず、システムが到達可能であることを確認します。つまり、全ての状態が何らかの入力によって到 達可能であることを保証します。 2. 変換: • 到達可能な状態空間表現を、特定の標準的な形式に変換します。この形式は通常、「Controllable Canonical Form」とも呼ばれます。 3. 特性の利用: • 到達可能正準形に変換することで、システムの特性がより直感的に理解しやすくなり、解析や制御設 計が容易になる場合があります。 到達可能正準形は、特に状態空間法を用いた制御システムの設計や解析において便利な形式への変換手法 として利用されます。

コメント(no.1) 58 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

おわり 59