概要
ついに，終着駅に向かってまいりました．
長かったですが，まだまださらにかかります・・・
H∞制御とは
いろいろなところで，語られているので軽く
min \|T_{zw}\|_{\infty}
を達成する，制御器$K$を設計し，制御に用いるということである．
計算方法：ガンマイタレーション
ガンマイタレーションという方法で計算する．LMIとかはまだ勉強していないが，有界実補題でLMIに変換すれば，可能だと思われる
次の定理を基本とした方法である．
\|T_{zw}\|_{\infty} <\gamma
を満たす制御器が存在すること
⇔次の3つが成立
(1) $H_{\infty} \in dom(Ric)$かつ，解$X_{\infty}$が正定
H_{\infty}=\left(
\begin{array}{ccc}
A &\gamma^{-2}B_1B_1^*-B_2B_2^*\\
-C_1^*C_1&-A^*\\
\end{array}
\right)
$\gamma = \infty$ならば，(1,2)成分が負定になり，可制御・可観測ならば，$dom(Ric)$になる(証明は余裕があるときに）．
(2) $J_{\infty} \in dom(Ric)$かつ，解$Y_{\infty}$が正定
J_{\infty}=\left(
\begin{array}{ccc}
A^* & \gamma^{-2}C_1^*C_1-C_2^*C_2\\
-B_1B_1^*&-A\\
\end{array}
\right)
(3)$\rho(A)$をAの最大ユークリッドノルムの固有値を示すとする．
\rho(A) = max_{1<i<n}\|\lambda_i\|
と定義する．この時，次を満たす．
\rho(X_{\infty}Y_{\infty} ) < \gamma^2
さらに，これらの条件が成立するとき，制御器の一つは
K_{sub}=\left(
\begin{array}{c|c}
\hat{A}_{\infty}&-Z_{\infty} L_{\infty}\\ \hline
F_{\infty}&0\\
\end{array}
\right)\\
\hat{A}_{\infty} = A+ \gamma^{-2}B_1 B_1^* X_{\infty} +B_2F_{\infty} + Z_{\infty}L_{\infty}C_2\\
F_{\infty} = -B_2^*X_{\infty}\\
L_{\infty} = -Y_{\infty}C_2^*\\
Z_{\infty} = (I-\gamma^{-2}Y_{\infty}X_{\infty})^{-1}\\
となる．
この定理を用いて，
min \|T_{zw}\|_{\infty}
を達成するには，
(1) 初期値$\gamma_1$を与え，3条件を満たすか確認する
(2) 満たす場合$\gamma_1$より小さい値で，3条件を満たすか確認する
(3) 満たさない場合，$\gamma_1$より大きい値で3条件を満たすか確認する
二分探索等の探索アルゴリズムを使用し，3条件を満たす，できるだけ小さい$\gamma$を計算する．
最後に，その時の$K_{sub}$を計算すれば，計算終了となる．ガンマ値を色々変えて繰り返すので，ガンマイタレーションと呼ばれる．
以降では，定理の証明をする．補題3つまず，証明し，定理を示す
ちなみに，すでに
\gamma > \| G \|_{\infty}
であると事と
R=\gamma ^2 I -D^*D\\
F=R^{-1}(B^*X+D^*C)\\
S=-C^*D\\
P=-C^*C\\
としたときのハミルトン行列H
H=
\left(
\begin{array}{ccc}
A+BR^{-1}S^*&-BR^{-1}B^*\\
-(P-SR^{-1}S^*)&-( A+BR^{-1}S^*)^*\\
\end{array}
\right)
が$H \in dom(Ric)$であることが
同値であることは示している．
Pythonで実装
証明は後にして，Scilabで設計したコントローラと，Pythonで自作した，ガンマイタレーションの結果の比較をしようと思います．
"""
二分探索で ガンマイタレーションを行い， 準最適制御器 を設計する．
"""
import numpy as np # Numpyライブラリ
import copy
def check(gamma,A,B1,B2,C1,C2):
"""
ハミルトン行列の生成
"""
H=copy.deepcopy(A)
hoge = np.dot(B1,B1.conj().T) / (gamma**2) - np.dot(B2,B2.conj().T)
H = np.append(H,hoge,axis = 1)
hoge2 = -1 * np.dot(C1.conj().T,C1)
hoge2 = np.append(hoge2,-1*A.conj().T ,axis =1 )
H = np.append(H,hoge2,axis = 0)
J=copy.deepcopy(A.conj().T)
hoge = np.dot(C1.conj().T,C1) / (gamma**2) - np.dot(C2.conj().T,C2)
J = np.append(J,hoge,axis = 1)
hoge2 = -1 * np.dot(B1,B1.conj().T)
hoge2 = np.append(hoge2,-1*A ,axis =1 )
J = np.append(J,hoge2,axis = 0)
"""
H
ハミルトン行列に固有値0を持たない，ハミルトン行列の解を求める，
"""
H=H.astype(np.float64)
Hw,Hv = np.linalg.eig(H)
for i in range(len(Hw)):
if Hw[i].real == 0:
print("H have pole 0")
return 0,0,0
v_num = 0
for i in range(len(Hw)):
if Hw[i].real < 0:
if v_num == 0:
V = Hv[:,i].reshape(H.shape[0],1)
v_num = v_num+1
else:
V = np.append(V,Hv[:,i].reshape(H.shape[0],1),axis = 1)
v_num = v_num+1
n=A.shape[0]
X1_H = V[0:n,0:n]
X2_H = V[2:n+2,0:n]
if np.linalg.det(X1_H) == 0:
print("H doesn't have Hojyuu seisitu")
return 1,0,0
else:
X_H = np.dot(X2_H,np.linalg.inv(X1_H))
"""
さらに，解の正定性を確認
"""
if np.all(np.linalg.eigvals(X_H) > 0):
print("X of H is positive")
else:
return 2,0,0
"""
リカッチ方程式の確認
"""
hoge = np.dot(A.conj().T,X_H)
hoge = hoge + np.dot(X_H,A)
R=np.dot(B1,B1.conj().T) / (gamma**2) - np.dot(B2,B2.conj().T)
hoge = hoge + np.dot(X_H,np.dot(R,X_H))
hoge = hoge + np.dot(C1.conj().T,C1)
#print(hoge)
#print(" is 0? Is riccati correct ?")
"""
J
ハミルトン行列に固有値0を持たない，ハミルトン行列の解を求める，
"""
J=J.astype(np.float64)
Jw,Jv = np.linalg.eig(J)
for i in range(len(Jw)):
if Jw[i].real == 0:
print("J have pole 0")
return 7,0,0
v_num = 0
for i in range(len(Jw)):
if Jw[i].real < 0:
if v_num == 0:
V = Jv[:,i].reshape(J.shape[0],1)
v_num = v_num+1
else:
V = np.append(V,Jv[:,i].reshape(J.shape[0],1),axis = 1)
v_num = v_num+1
n=A.shape[0]
X1_J = V[0:n,0:n]
X2_J = V[2:n+2,0:n]
if np.linalg.det(X1_J) == 0:
print("J doesn't have Hojyuu seisitu")
return 3,0,0
else:
X_J = np.dot(X2_J,np.linalg.inv(X1_J))
"""
さらに，解の正定性を確認
"""
if np.all(np.linalg.eigvals(X_J) > 0):
print("X of J is positive")
else:
return 4,0,0
"""
リカッチ方程式の確認
"""
hoge = np.dot(A,X_J)
hoge = hoge + np.dot(X_J,A.conj().T)
R=np.dot(C1.conj().T,C1) / (gamma**2) - np.dot(C2.conj().T,C2)
hoge = hoge + np.dot(X_J,np.dot(R,X_J))
hoge = hoge + np.dot(B1,B1.conj().T)
#print(hoge)
#print(" is 0? Is riccati correct ?")
"""
スペクトル半径が gamma^2 以下である事の確認
"""
hoge = np.dot(X_H,X_J)
spectro_radius = max(np.abs(np.linalg.eigvals(hoge)))
if spectro_radius < gamma**2:
print("spectro radius < gamma^2")
return 6,X_H,X_J
else:
print("spectro radius error")
return 5,0,0
A=np.array([[-1,10],[2,-0.1]])
w,v = np.linalg.eig(A)
for i in range(len(w)):
if w[i].real < 0:
print("A is unstable")
if w[i].real == 0:
print("A have pole 0")
B1=np.array([[1,0],[2,0]])
B2=np.array([[1],[2]])
C1=np.array([[1,3],[0,0]])
C2=np.array([[1,3]])
D12 = np.array([[0],[1]])
D21= np.array([[0,1]])
"""
仮定の確認
"""
hoge = np.dot(D12.conj().T,D12)
I=np.diag(np.ones( np.shape(D12)[1] ))
if ( hoge == I).all:
print("D12* D12 = I OK")
else:
print("D12* D12 = I error")
hoge = np.dot(D21,D21.conj().T)
I=np.diag(np.ones( np.shape(D21)[0] ))
if (hoge == I).all():
print("D21 D21* = I OK")
else:
print("D21 D21* = I error")
hoge = np.dot(D12.conj().T,C1)
zeros = np.zeros(hoge.shape)
if (hoge == zeros).all():
print("D12* C1 = 0 OK")
hoge = np.dot(B1,D21.conj().T)
zeros = np.zeros(hoge.shape)
if (hoge == zeros).all():
print("B1 D21* = 0 OK")
"""
可制御・可観測の確認
"""
VC = copy.deepcopy(B1)
hoge = copy.deepcopy(A)
for i in range(A.shape[0]-1):
hoge = np.dot(hoge,B1)
VC=np.append(VC,hoge,axis = 1) #横につなげる
if A.shape[0] == np.linalg.matrix_rank(VC):
print("A B1 is controllable")
VC = copy.deepcopy(B2)
hoge = copy.deepcopy(A)
for i in range(A.shape[0]-1):
hoge = np.dot(hoge,B2)
VC=np.append(VC,hoge,axis = 1) #横につなげる
if A.shape[0] == np.linalg.matrix_rank(VC):
print("A B2 is controllable")
VO = copy.deepcopy(C1)
hoge = copy.deepcopy(A)
for i in range(A.shape[0]-1):
hoge = np.dot(C1,hoge)
VO=np.append(VO,hoge,axis = 0) #横につなげる
if A.shape[0] == np.linalg.matrix_rank(VO):
print("C1 A is observable")
VO = copy.deepcopy(C2)
hoge = copy.deepcopy(A)
for i in range(A.shape[0]-1):
hoge = np.dot(C2,hoge)
VO=np.append(VO,hoge,axis = 0) #横につなげる
if A.shape[0] == np.linalg.matrix_rank(VO):
print("C2 A is observable")
"""
ガンマイタレーションの開始
"""
gamma_up = 10**10
gamma_down = 0
while True:
if gamma_up - gamma_down <0.001:
break
gamma = (gamma_up - gamma_down)/2 + gamma_down
flag,X_H,X_J = check(gamma,A,B1,B2,C1,C2)
if flag == 6: #条件を満たす
gamma_up = gamma
else: #条件を満たさない
gamma_down = gamma
print(gamma,flag )
print("gamma_up",gamma_up)
print("gamma_down",gamma_down)
gamma = (gamma_up - gamma_down)/2 + gamma_down
F = -1*np.dot(B2.conj().T,X_H)
n=A.shape[0]
eye = np.diag(np.ones(n))
Z = np.linalg.inv(eye - np.dot(X_J,X_H)/(gamma**2))
L = -1*np.dot(X_J,C2.conj().T)
A_K = A + np.dot(B1,np.dot(B1.conj().T,X_H))/(gamma**2)+np.dot(B2,F)+np.dot(Z,np.dot(L,C2))
B_K = -1*np.dot(Z,L)
C_K = F
D_K = 0
print("A_K",A_K)
print("B_K",B_K)
print("C_K",C_K)
print("D_K",D_K)
scilab
A=[-1 10
2 -0.1]
B=[1 0 1
2 0 2]
C=[1 3
0 0
1 3]
D=[0 0 0
0 0 1
0 1 0]
x0=[0
0.05]
sl = syslin('c',A,B,C,D,x0);
gopt = gamitg(sl,[1,1])
//[K,ro] = h_inf(sl,[1,1],0,100,100)
K=ccontrg(sl,[1,1],100);
ガンマの最適値はほぼ一致しますが，設計されるコントローラが一致しません・・・
やめます
途中ですが，こんなことを頑張る意味がない気がしてきました．
理論がわからなくても何とかなりますし，時間をかける必要がない気がしてきました．
もう，やめます．