Notebook

微分代数方程式系 (DAEs) の動的最適制御を行う¶

JModelica.org Users Guide - Version 2.10 6.2.A first example に掲載されている微分代数方程式の動的最適制御問題 (dynamimc optimal control problem) を実習します。

概要¶

制御対象の微分代数方程式モデル (DAEs)¶

Van der Pol oscillator のモデルです。状態変数を $x_1$ , $x_2$ , 制御変数を $u$ とします。状態変数の初期値を設定すると、以下の微分代数方程式系に従って状態変数が変化します。

$\begin{align} \frac{d x_1}{dt} &= (1-x_2^2)x_1-x_2+u, \\ \frac{d x_2}{dt} &= x_1 \end{align}$

動的最適制御問題 (dynamic optimal control problem)¶

制約条件 (constraint) $u \leq 0.75$ のもとで制御変数を調整し、 startTime から finalTime まで状態変数を移動させたときに、目標関数

$\int_{startTime}^{finalTime}{x_1^2 + x_2^2 + u^2}dt$ が最小になるようにします。

モデルの作成¶

最適化モデル VDP_Opt のソースコードを示します。 objectiveIntegrand 目標関数の被積分関数を設定します。

In [1]:

%%writefile VDP_Opt.mop
optimization VDP_Opt (objectiveIntegrand = x1^2 + x2^2 + u^2,
                      startTime = 0,
                      finalTime = 20)
  // The states
  Real x1(start = 0, fixed = true);
  Real x2(start = 1, fixed = true);
  // The control signal
  input Real u;
equation
  der(x1) = (1 - x2^2) * x1 - x2 + u;
  der(x2) = x1;
constraint
  u <= 0.75;
end VDP_Opt;

Overwriting VDP_Opt.mop

最適化モデルをロードします。

In [2]:

from pyjmi import transfer_optimization_problem
op = transfer_optimization_problem("VDP_Opt", "VDP_Opt.mop")

動的最適制御シミュレーションを行います。

In [3]:

res = op.optimize()

Initialization time: 0.20 seconds

Total time: 0.25 seconds
Pre-processing time: 0.00 seconds
Solution time: 0.03 seconds
Post-processing time: 0.02 seconds

シミュレーション結果を取得してプロットします。

In [4]:

x1 = res['x1']
x2 = res['x2']
u = res['u']
t = res['time']

In [5]:

%matplotlib notebook
import matplotlib.pyplot as plt
import numpy as np
plt.figure(1)

plt.subplot(311)
plt.plot(t,x1)
plt.xlim(0,20)
plt.ylim(-0.5,0.1)
plt.xticks(np.arange(0,20,step=5),color="None")
plt.yticks(np.arange(-0.5,0.2,step=0.1))
plt.grid()
plt.ylabel('$x_1$')

plt.subplot(312)
plt.plot(t,x2)
plt.xlim(0,20)
plt.ylim(-0.2,1.0)
plt.xticks(np.arange(0,20,step=5),color="None")
plt.yticks(np.arange(-0.2,1.2,step=0.2))
plt.grid()
plt.ylabel('$x_2$')

plt.subplot(313)
plt.plot(t,u)
plt.xlim(0,20)
plt.ylim(-0.6,0.8)
plt.xticks(np.arange(0,20,step=5),color="None")
plt.yticks(np.arange(-0.6,1.0,step=0.2))
plt.grid()
plt.ylabel('$u$')

Out[5]:

Text(0,0.5,u'$u$')

In [ ]: