#!/usr/bin/env python
# coding: utf-8

# ## 微分代数方程式系 (DAEs) の動的最適制御を行う
# [JModelica.org Users Guide - Version 2.10](https://jmodelica.org/downloads/UsersGuide-2.10.pdf) 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]:


get_ipython().run_cell_magic('writefile', 'VDP_Opt.mop', 'optimization VDP_Opt (objectiveIntegrand = x1^2 + x2^2 + u^2,\n                      startTime = 0,\n                      finalTime = 20)\n  // The states\n  Real x1(start = 0, fixed = true);\n  Real x2(start = 1, fixed = true);\n  // The control signal\n  input Real u;\nequation\n  der(x1) = (1 - x2^2) * x1 - x2 + u;\n  der(x2) = x1;\nconstraint\n  u <= 0.75;\nend VDP_Opt;\n')


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

# In[2]:


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


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

# In[3]:


res = op.optimize()


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

# In[4]:


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


# In[5]:


get_ipython().run_line_magic('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$')


# In[ ]: