平田氏の著書「実践ロバスト制御」↓を無料のPython-Controlで勉強する試み⑦。
今まで，教科書のMATLABサンプルの通り，hinfsynを使ってH∞制御問題を解いてきたが，mixsynというコマンドもあって，使い方を調べてみたら，混合感度問題を解くならこっちを使った方が楽だった。
前の記事↓で書いたように，
manva.hatenablog.com
hinfsynを使う場合のH∞制御系の設計手順は以下の通り。
hinfsynを使う場合の設計手順
１．制御対象の摂動を定義する
２．乗法的摂動のボードゲイン線図を描く
３．それを覆うように重みを決める
４． 感度関数に対する重みを決める
５．一般化プラントGを定義
６． 解く！ hinfsyn
mixsynを使う場合，上記の5番がいらない。
mixsynを使う場合の設計手順
１．制御対象の摂動を定義する
２．乗法的摂動のボードゲイン線図を描く
３．それを覆うように重みを決める
４． 感度関数に対する重みを決める
５．解く！ mixsyn
Python-Controlのcontrol.mixsynへの引数は(g,w1,w2,w3)で，gは一般化プラントではなくてノミナルの制御対象を与え，w1に感度関数に対する重み，w3に相補感度関数に対する重み
を与える。
control.mixsyn — Python Control Systems Library dev documentation
基本はMATLABの仕様に合わせているので，使い方はMATLABのドキュメント↓を参照するとよい。
Mixed-sensitivity H∞ synthesis method for robust control loop-shaping design - MATLAB mixsyn - MathWorks 日本
↓の図がわかりやすい。
Mixed-Sensitivity Loop Shaping - MATLAB & Simulink - MathWorks 日本
w2は今回はいらないのだが，
K, CL, info = control.mixsyn(Pn,Ws,None,Wt)
としたら応答が返ってこなかったので，
K, CL, info = control.mixsyn(Pn,Ws,tf(0.001,1),Wt)
のように伝達関数で小さい値を入れておいたらできた。
↓の問題をmixsynで解くと
無料で実践ロバスト制御④ ーH∞制御問題を解くー - manvaのエンジニアリング魂
↓のような感じ。
# 実践ロバスト制御4.2節 import numpy as np import matplotlib.pyplot as plt from control.matlab import * import control # %% パラメータの定義 m1 = 0.8; m2 = 0.2; k1 = 100; k2 = 300; c1 = 1; c2 = 0.3; Ks = 100; # %% 運動方程式のM,K,C行列を定義 M = np.matrix([ [ m1, 0 ], [ 0, m2 ] ]); C = np.matrix([ [ c1+c2, -c2 ], [ -c2, c2 ] ]); F = np.matrix([ [ Ks ], [ 0 ] ]); K = np.matrix([ [ k1+k2, -k2 ], [ -k2, k2 ] ]); # %% 状態空間実現 iM = np.linalg.inv(M); Ap = np.concatenate([ np.concatenate([np.zeros((2,2)), np.eye(2)], axis=1), np.concatenate([ -iM*K, -iM*C ], axis=1) ]); Bp = np.concatenate([ np.zeros((2,1)), iM*F ]); Cp = [ 0, 1, 0, 0 ]; Dp = 0; # %% ノミナル制御対象の定義 P_nominal = ss( Ap, Bp, Cp, Dp ); Pn = ss2tf(P_nominal); #%% 重み関数の定義 s = tf('s'); Wm = 3*s**2/(s**2+18*s+45**2); Ws = 15/(s + 1.5e-2); #% Ws Wt = Wm; #% Wt #%% H-inf制御器の計算（実行4.1） K, CL, info = control.mixsyn(Pn,Ws,tf(0.001,1),Wt) #混合感度問題 print( info ) plt.figure(4) mag1, phase1, om1 = bode( P_nominal, plot=False ); plt.semilogx( om1, 20*np.log10(mag1), '--', label='$P_{nominal}$' ); mag2, phase2, om2 = bode( K, plot=False ); plt.semilogx( om2, 20*np.log10(mag2), label='$K$' ); plt.ylim(-60, 60); plt.xlim(1e0, 1e2); plt.legend()
実行結果
mixsynでも教科書の図4.8に相当する結果が得られた。
修正混合感度問題の場合
修正混合感度問題の場合は，感度関数をとおけばよい。
↓の問題をmixsynで解くには，
無料で実践ロバスト制御⑤ ー修正混合感度問題ー - manvaのエンジニアリング魂
上記プログラムでmixsynのところを↓のように書き換えるだけでよい。
WpsP = Ws*0.8*Pn; K, CL, info = control.mixsyn(Pn,WpsP,tf(0.001,1),Wt)
実行結果
教科書の図4.14に相当する結果が得られた。
修正混合感度問題 h2synの場合
ついでに，どんなものかよくわかっていないが，h2synも試してみた。こちらはhinfsynと同様の使い方で，一般化プラントを与える。↓の記事で
無料で実践ロバスト制御⑤ ー修正混合感度問題ー - manvaのエンジニアリング魂
移植した「プログラム4.5」のhinfsynのところを，↓のように書き換えるだけでできる。
K = control.h2syn(G,1,1)
