11回目（2020年度秋）　応用プログラミング演習（制御工学）(1)

制御シミュレーションの数値解法

これから学ぶであろう，制御工学について情報処理演習の観点から扱ってみたいと思います．
制御工学で学ぶのは，線形系であるため，基本は解析解が求められる状態を扱います．(紙と鉛筆で解ける)
しかしながら，昨今の制御系はディジタルディバイスを利用するため，連続量でなくセンサなどから得られる「離散量を基にした制御」を考えておく必要があります．

制御理論のごく基礎(知っていて損はありません)

私たちの身の周りには，空調の温度設定や，自動車の自動速度設定など，機器などの動作状態を指示通りにコントロール（＝制御）したい場合が多く存在します．
何かの対象の状態を計測し，それを元に所望の設定通りに状態を制御する方法の基本は，フィードバック制御と呼ばれる方法です．
動作原理はいたってシンプルで，

制御対象の「状態」をセンサなどで計測
「設定値」と計測値を比較して偏差（ズレ）を算出
偏差をもとに，対象に対して何らかの作用(修正)

以下に原理図を示します．

上記の図の記号は，

SV(Set Variable)：目標値，設定した値（エアコンの設定温度など）
PV(Process Variable)：計測値，センサの出力値（温度センサの出力など）
MV(Manupilative Variable)：操作量，制御器（コントローラ）が作り出す指令（ヒートポンプの出力の増減など）

となっており，目標値と計測値を比較し，その差（偏差）を解消するように制御を行うのが基本の機能となっています．
制御対象の特性に合わせて制御器を工夫することにより制御特性が改善されます．

操作に対して，対象が即座に反応して設定値に到達するならば，制御の必要はないが，現実には操作と反応には時間差が生じます．
例えば，大きな部屋の空調温度変化（熱容量）や，回転モーメントの大きな回転運動系（慣性）などの物理現象がこれにあたるのです．

モータ回転数の制御

制御対象として一般的に用いられるモータを考えましょう．
モータの応答はかけた電圧に対して「1次遅れ」と呼ばれる応答を示すことが知られています．
以下に1次遅れの例を図示し，スイッチを入れたような入力に対する応答をステップ応答と呼んでいます．

これを数式で示すと，

y(t) = K (1 - exp(-t / T))

となるものであり，Kはゲイン係数（定数），Tは時定数と呼ばれるもので，1次遅れの基本的な特性を示しています．時定数Tは応答の遅れ具合と考えればよいのです．
ある値に向かってはじめは急に応答し，徐々に穏やかに回転が上昇します．身近なものでは，掃除機のモータ音などを思い出してみれば分かりやすいかと思います．

ディジタル制御

制御理論を学ぶときにはアナログ量を基にした解析を行いますが，現在はマイコンを主体とした制御装置を利用することが大半ですので，実際にはディジタル制御を行います．
ディジタル制御を行うため，周期的な演算を繰り返します．制御周期Δｔを決め，その間に，計測，演算，操作量MVの出力を周期的に繰り返して実行します．
制御システムの機能は，設定値 SV および計測値 PV から，MVをどのように決定するかが問題となります．

従って，Δtは短い時間とはいえ，間隔をおいた断続的な制御を行っているのです．
アナログ量では時間の関数であり，PV(t), MV(t), e(t)と示すことができましたが，ディジタルの場合，Δtごとにn周期目と考えることができるので，PVｎ，MVｎ，enのように考えます．

MVｎ=MVｎ-1+ΔMVｎ

この式は，前回の操作量 MVn-1 に今回の変化分 ΔMVn だけ加える，という意味である．
1周期前でしたらPVｎ-1，MVｎ-1，en-1となるわけです．

比例(P)制御コントローラ
「制御器」の仕組みを調整することで，システム全体の応答を様々に変化させることが可能である．
最も基本的な制御を，比例(Propotional)制御，または単にP制御と呼び，制御器の感度を調整することで，応答の様子を変化させることが可能です．
P制御の他にも，PD(Proportional-Differential)制御，PID(Proportional Integral Differential)制御などがよく用いられます．
偏差e(入力-応答)に比例ゲインKpをかけて操作量MVを生成します．これは，

MV = Kp e = Kp (y(t) -SV)

で示されます．Kpを調整すれば，偏差に対して敏感に修正する特性にすることも可能となります．
偏差に対して，ばねのように誤差を減らそうと機能します．Kpがちょうどばね定数のようなものとして考えて下さい．
先の表記に則り，検出器が応答y(t)をそのまま示すとすれば，

MV = Kp e = Kp (PV -SV)

と示すことができます．

位置型ディジタル制御アルゴリズム

さて，アナログ量で扱う場合での時刻tにおける偏差と操作量MV(t)はそれぞれ，

e(t) = SV(t) - PV(t)，　MV(t) = Kp e(t) = Kp (SV(t) - PV(t))

となりますが，ディジタル制御では，現在の周期nで考えますので，

en = SVn - PVn，　MVn = Kp en + MV0 = Kp (SVn - PVn) + MV0

と考えればよいのです．制御を開始した時の初期値MV0を加える必要があります．
すなわち，P制御においては，ΔMVnは偏差の変化量と考えることができるので，

ΔMVn = en - en-1

と計算します．

速度型ディジタル制御アルゴリズム

同様に詳細は省略しますが，計算量を少なくする制御アルゴリズムに速度型制御アルゴリズムと呼ばれるものがあり，以下の式で操作量を示すことが可能です．

MVn = MVn-1 + ΔMVn --- (a)

この式の意味は，前回の操作量MVn-1に今回の変化分ΔMVnだけを加える，という意味です．

P制御においては，偏差の変化量と考えることができるので，

ΔMVn = Kp (en - en-1) --- (b)

と計算することができます．

※　P制御では位置型でも速度型でも計算量に大きな違いはありませんが，他の制御方式で，積分が入ったり，フィルタが組み合わされるような場合に効果があることもあります．

1次遅れモデルの離散表現

アナログ値では1次遅れ系においてTは時定数で，目標値と操作量は以下の関係にあります．

(d/dt)PV + (1/T) PV = MV

周期Δｔごとに制御するため，n周期後の応答は以下の式で示されます．

PVn = (T / (T + Δｔ)) PVn-1 + (Δｔ / (T + Δｔ)) MVn --- (c)

なお，P制御をかけているので，今与えるべき操作量は，

MVn = Kp en = Kp (PVn-1 - SVｎ) --- (d)

参考）位置型，速度型におけるMVnの与え方

位置型であれば，式(d)を式(c)に代入すれば，現在の制御周期における応答が計算できます．
速度型に関しては，前回の制御周期における偏差との差を計算する式(b)を計算し，前回の制御周期における操作量に加える式(a)を行った上で式(c)に代入すれば同様に計算できます．

速度型の方がひと手間多いように見えますが，前回の制御周期との差だけで計算できるので，初期状態からのデータの蓄積に依存しないで済み，その結果，制御法を途中で切り替えた時に急に値が変化することなくスムースに移行できるなどのメリットがあります．
あるいは離散化の際には必ず誤差が発生しますが，誤差の蓄積を抑えられる可能性もある点で利用されています．

シミュレーションプログラム作成にあたって

上記のように速度型制御アルゴリズムを扱っているため，以下のような変数のみを用意すれば
制御できます．

目標値SV, 制御対象の時定数T, 制御周期Δt,　比例ゲインKp
現制御周期の偏差eｎ, 前制御周期の偏差eｎ-1
現制御周期の操作量MVｎと変化分ΔMVn, 前制御周期の操作量MVn-1
現制御周期の応答PVn

を用意すればよいのです．nやn-1と言っても全周期分用意する必要はなく，現制御周期と前制御周期に関するデータを一時的に記録しておけばよく，fprintfを用いてファイルに書き出してしまえば良いのです．

もちろん配列に記録してからファイルに書き出してもかまいません．

以上が，制御に関するシミュレーションについての説明です．かなり簡素化していますが，プログラムもやってみるとそれほど大変ではないことが分かると思います．

（参考）　関数ポインタ
変数を関数で渡す方法は既に学んでいると思いますが，実は関数もポインタを利用してやりとりをする
ことが可能なのです．
以下に例を示したいと思います．

double fnc(double x){
　return x;
}

main(){
　double x;
　double (*pfnc)(double);
　pfnc=fnc;

　x=10.0;

　*pfnc(x);
}

通常の関数を宣言します．

メイン関数内で関数ポインタを宣言します．
先の関数のアドレスを関数ポインタに送ります．

関数ポインタに引数を渡し，先の関数を実行します．

または，

double fnc(double x){
　return x;
}

main(){
　double x;
　double (*pf)(double);
　*pf=fnc;

　x=10.0;

　
}

通常の関数を宣言します．

メイン関数内で関数のポインタを宣言します．
先の関数を関数のポインタの中身に送ります．

のようにしてもやりとりができます．
関数をわざわざポインタにする理由は何かというと，

関数の再利用性を上げることができる
関数を動的に何度でも呼び出して利用できる
関数に対して異なった種類の関数を呼び込むことができる．

既に作っておいた関数を利用する際，ポインタで受け渡しをすることで，例えば関数名が合わなくてもその制約をなくして利用できるので再利用性が向上します．また，関数をモジュール（ユニット）のようなものとして利用できるので関数内のアルゴリズムの構造を何度も形を変えて利用できるという点があります．何かの数値を算出するための関数ではなく，計算するための仕組みとしての関数にでき，汎用性，再利用性が格段に上がります．

授業の課題であれば，プログラムが小規模な上に，あまり再利用する機会がないのでメリットがつかみにくいと思いますが，自作のプログラムを作る際は大変便利ですので参考にしてみて下さい．

関数ポインタ（補足）
関数ポインタを関数に呼び込む方法が分かりにくいと思いますので，追加説明します．以下の例を見て下さい．
fnc_1の引数として関数ポインタを用いています．メイン関数内では関数fnc_A, fnc_Bをこのfnc_1に渡しています．
fnc_1の引数(*fnc)(double)は関数ポインタの中身を表しているので，(*fnc) = fnc_Aとして関数ごと引き渡すことができると考えて下さい．

#include <stdio.h>

double fnc_1(double (*fnc)(double), double a){
　return *fnc(a);
}

double fnc_A(double a){
　return a;
}

double fnc_B(double a){
　return a*a;
}

main(){
　double h=2;

　printf("fnc_A=%f\n", fnc_1(fnc_A, h));
　printf("fnc_B=%f\n", fnc_1(fnc_B, h));
}

(参考)：　時間ごとの制御処理（タイマ割り込み）

如何でしょうか？今回作成した制御プログラムはシミュレーションのためのものであり，データをファイルに出力したり，画面に表示するだけかと思うかも知れませんが，センサからのテータを取得し，モータに命令を送るための信号生成とそれを受ける回路を用意すれば，今回の授業で扱ったアルゴリズムを使うことで実際の制御を行うことができるのです．

プログラムにすると以下のような流れとなります．

int main(){
　タイマ初期化関数
　タイマ（カウンタ）スタート関数

　while(1){
　　if(変数==○){
　　　割り込み用処理
　　　カウンタ変数=0;
　　}
　通常の処理
　}
}

int タイマ割り込み本体関数(){
　タイマ停止
　割り込み処理
　カウンタリセット
　タイマスタート
}

main関数
　タイマ初期化
　タイマ（カウンタ）スタート

　ループ
　　割り込み発生（カウンタ変数を満たす）
　　　割り込み用処理の実行
　　　カウンタ変数のリセット

　割り込み以外の通常処理

タイマ割り込み本体関数
　タイマ停止
　割り込み処理
　カウンタリセット
　タイマスタート

これは割り込みという制御で，特に機械類の制御の場合，短い時間の周期（1～50msec）位の制御周期を保ちながら動作する必要があるのです．

これまでのプログラムのように，一つの処理が終わったら次の処理，但し，それぞれ完了する時間は分からない，という状態では制御ができないので，タイマークロックを基準にし，時間が来たら必ず動作させるような処理を行います．

これを図で示すと以下のようになります．

課題1

上記の位置型アルゴリズムを基に，速度法の関数を作成して，制御シミュレーションを行う．

関数1：MVの計算

式(d)より，現時刻の制御周期における操作量MVnを生成する．

MVn = Kp en = Kp ( SV - PVn-1)　　 --- (d)

関数の計算に必要な値を引数として、MVを計算して返す関数を作成する．

関数2：PVの計算

式(c)より，時定数 T のシステムの現在の応答PVを計算する．
これは制御周期 Δt の時系列データとなるので，ループ処理とする． PVnを画面に順次表示するか，またはファイルに書き出す．

この関数はシミュレーション終了時刻 Te を受け取る必要がある．

PVn = ( T / ( T + Δt )) PVn-1 + ( Δt / ( T + Δt )) MVn　　--- (c)

ひとつ前の周期の値 PVn-1 が必要であることに注意．

課題2

時定数T=2.0秒，シミュレーション時間SimTは5～10.0秒にし，目標値はSV=5000[rpm]，Kpを自由に調整して良いと思われる応答になるよう調整し,Excelでグラフを作ろう．

制御周期 Δt は2-3種類，Kp も2-3種類変え，比較のグラフを作成すること．どのような相違が生ずるか考察せよ．
（すべてのグラフを重ねると見づらいので， Δt ごとに，Kp違いの PV 値を1つのグラフ内に描くと良い．）
なるべく速く目標値に到達するのが，制御の目標であるので，Kを調整してグラフの変化を比較しよう．

処理の流れ

１シミュレーションの各定数を決める
２ループを設定し，刻み時間 Δt 間隔で
-操作量 MVnの値を計算する．
-応答 PVnの値を求め，ファイルや画面に書き出す，または配列に格納．
３シミュレーション終了時刻までループを繰り返す．

により計算できる．なお，初期値PV0以前は 0 と置き，前回の応答値PVn-1 は別の変数に一時的に格納すればよい．

Kp を0.5， 20， 50 など，大幅に変えて比較すると良い．，

注：比例制御だけでは偏差（オフセット）が0にならず，制御量を目標値に一致させることができないことが理論的に知られている．これを定常偏差と言う．

参考ソースファイル：ctrl_sample.cpp

注：参考ソースファイルは理解を促すためのものであり，丸写しでの提出は厳禁．

ヒント（以前の課題ですが，ソースコードが参考になれば使ってみて下さい．）
比例制御について以下のファイルをダウンロードし，シミュレーションを行いましょう．どちらを用いても構いません．

速度型アルゴリズム　p_p_cl.c
位置型アルゴリズム　p_v_cl.c

時定数T=2.0秒，シミュレーション時間SimTは5～10.0秒にし，目標値はSV=5000[rpm]，Kpを自由に調整して良いと思われる応答になるよう調整し．グラフに出力しましょう．

いくつか比較できればなおよいので挑戦してみて下さい．
ソースファイルは完全ではありません．このファイルは開発途中の様相となっているため，余計な機能があったり，書きかけの記述があったり，コメントアウトがあったりするので，適宜整えて必要なソースを作成して下さい．画面に数値を出力したり，ファイルに数値が書き出せるよう，変更して下さい．
なお，エクセルに読み込ませるには，以下のような形式で記述するとよいので参考にして下さい．

エクセルに読み込めるファイルの形式の代表(これらをまとめてCSVファイルともいう)

ファイル出力例
CSV（カンマ区切り）ファイル	ＴSV（タブ区切り）ファイル	SSV（スペース区切り）ファイル
1.00000, 2.00000, 3.00000	1.00000 2.00000 3.00000	1.00000 2.00000 3.00000
プログラムの表記例
printf("%lf,%lf,%lf\n", a, b, c);	printf("%lf\t%lf\t%lf\t\n", a, b, c);	printf("%lf %lf %lf\n", a, b, c);

本サンプルファイルは，レトロなコマンドライン上でグラフ概形を表示する機能を入れています．

提出するもの：　自ら書き換えたソースファイル，Excelなどで作成したグラフの入ったファイル

参考
PID制御について以下のファイルをダウンロードし，シミュレーションを行いましょう．どちらを用いても構いません．

速度型アルゴリズム　pid_p_cl.c
位置型アルゴリズム　pid_v_cl.c

時定数T=2.0秒，シミュレーション時間SimTは5～10.0秒にし，Kpを自由に調整して良いと思われる応答になるよう調整し．グラフに出力しましょう．
ソースファイルは完全ではありません．画面に数値を出力したり，ファイルに数値が書き出せるよう，変更して下さい．

提出するもの：　自ら書き換えたソースファイル，Excelなどで作成したグラフの入ったファイル

参考　制御を理解するためのチャレンジ課題

上記のファイルはコマンドライン上でグラフ概形を知るための苦肉の策であるため，本当に概形でしかないので，グラフィックを利用することが有効です．以下のページを読んでおき，もし可能であれば，本日の課題を改変して，グラフィック表示できるようにしてみてください．

参考ページ　windows グラフィック　入門

さらに，もし可能であれば，引数に関数ポインタを使ってみてください．

関数ポインタについて