# ソルバを使う

**MatrixEquations.jl**を使って最適レギュレータを設計する．  
公式ページ : [https://docs.juliahub.com/MatrixEquations/1uOBF/1.1.4/](https://docs.juliahub.com/MatrixEquations/1uOBF/1.1.4/)


```julia
using Plots
using LinearAlgebra
using MatrixEquations
```

システムを定義  

$$
    \begin{align}
    \dot{x} &= Ax+Bu\\
    y &= Cx
    \end{align}\\
    A\in\mathbb{R}^{2\times2},~B\in\mathbb{R}^{2\times2},~C\in\mathbb{R}^{2\times2}
$$


```julia
x₀ = [1, 0.5]
A = [
    1.1 2
    -0.3 -1
]
B = [
    1 2
    0.847 3
]
C = [
    1. 0.
    0. 1.
];
```

安定性を調べてみる


```julia
eigvals(A)
```


    2-element Vector{Float64}:
     -0.6588723439378912
      0.7588723439378913


負の固有値を持つので不安定．  
初期値応答を調べる．  


```julia
t = 0:0.01:5
x = Vector{typeof(x₀)}(undef, length(t))
for i in 1:length(t)
    x[i] = exp(A .* t[i]) * x₀
end

x = permutedims(hcat(x...))
plot(t, x[:, 1], label="x1")
plot!(t, x[:, 2], label="x2")
plot!(xlabel="t")
savefig("lqr.png")
```

  
発散することがわかる．  

## リカッチ方程式を解く
`MatrixEquations`の`arec`を使う．  
ドキュメント : [https://docs.juliahub.com/MatrixEquations/1uOBF/1.1.4/autodocs/#MatrixEquations.arec](https://docs.juliahub.com/MatrixEquations/1uOBF/1.1.4/autodocs/#MatrixEquations.arec)  

`arec`は次の連続代数リカッチ方程式を解く．  




```julia
Q = diagm([10.0, 10.0])
R = diagm([1., 1.])
S = zero(B)
P, E, F = arec(A, B, R, Q, S);
```

戻り値は次の通り．  
* `P` : リカッチ方程式の解  
* `E` : 最適レギュレータの閉ループの固有値  
* `F` : 最適フィードバック行列  


```julia
P  # リカッチ方程式の解
```


    2×2 Matrix{Float64}:
      3.34568  -1.07266
     -1.07266   1.30235



```julia
E  # 閉ループの固有値
```


    2-element Vector{Float64}:
     -11.862446758953242
      -2.732479989130699



```julia
F  # 最適フィードバック行列
```


    2×2 Matrix{Float64}:
     2.43714  0.0304348
     3.47339  1.76174


最適フィードバック行列を入力に使ったときのシステムの時間発展を調べる．  


```julia
Ā = A .- (B * F)  # 新しいA行列

x2 = Vector{typeof(x₀)}(undef, length(t))
u2 = Vector{typeof(x₀)}(undef, length(t))
for i in 1:length(t)
    x2[i] = exp(Ā .* t[i]) * x₀
    u2[i] = -F * x2[i]
end

x2 = permutedims(hcat(x2...))
figx = plot(t, x2[:, 1], label="x1")
plot!(figx, t, x2[:, 2], label="x2")
plot!(figx, xlabel="t")
savefig(figx, "lqr2.png")

u2 = permutedims(hcat(u2...))
figu = plot(t, u2[:, 1], label="u1")
plot!(figu, t, u2[:, 2], label="u2")
plot!(figu, xlabel="t")
savefig(figu, "lqr2u.png")
```

  
  
  
安定化されていることがわかる．  

3810abe5f68090d482f9ab32a33d878bb8679dbf
ipynb
Jupyter Notebook
lqr/julia_src/using_solver.ipynb
YoshimitsuMatsutaIe/abc_2022
9c6fb487c7ec22fdc57cc1eb0abec4c9786ad995
MIT
lqr/julia_src/using_solver.ipynb
YoshimitsuMatsutaIe/abc_2022
9c6fb487c7ec22fdc57cc1eb0abec4c9786ad995
MIT
lqr/julia_src/using_solver.ipynb
YoshimitsuMatsutaIe/abc_2022
9c6fb487c7ec22fdc57cc1eb0abec4c9786ad995
MIT
Qwen/Qwen-72B
1. YES
2. YES