[PYTHON] Solve the initial value problem of ordinary differential equations with JModelica

Purpose / aim

• I want to use the trendy 1DCAE
• I want to check the usability of the interface of JModelica
• Study the Modelica language
• I want to find the good points of JModelica before being told that OpenModelica is okay

For the environment construction method, refer to the article Installing JModelica on Ubuntu.

Basic usage of JModelica

1. Prepare a Model file (* .mo)
2. Compile to Function Mockup Units (FMUs)
3. Load FMUs
4. Calculate the loaded FMUs
5. View results

Of course, it is also possible to read compiled FMUs from step 3.

Solve ordinary differential equations

According to Google Sensei, when you say Hello World in Modelica, there are many initial value problems of linear ordinary differential equations of the first order.

\begin{eqnarray}
\frac{dx(t)}{dt} &=& -x(t) \\
  x(0) & = & 1
\end{eqnarray}

The method of solving with paper and pencil is left to the textbook, and the derivation of the analytical solution is done with SymPy. Execute the following under the environment where SymPy can be used.

python

import sympy
x = sympy.Function("x"); t,C1 = sympy.symbols("t C1")
#x(t)Solve about x(t) == C1*exp(-t)
ans = sympy.dsolve(x(t).diff(t)+x(t),x(t))
#Calculate the constant of integration C1(t=0,x(0)=1)And substitute for the expression of ans
C = {C1:ans.subs(x(t),1).subs(t,0).lhs}
ans.subs(C)
#--> x(t) == exp(-t)

From the above, the analytical solution is as follows.

\begin{eqnarray}
  x(t) & = & \exp(-t)
\end{eqnarray}

Solve ordinary differential equations with JModelica

1. Preparation of model file

Prepare the following model file

ode_test.mo

model HelloWorld 
	Real x(start=1);
equation
	der(x)= -x;
end HelloWorld;

Lines 1-5: Model (class) definition 2nd line: Definition of state variable x with initial value 1 3rd line: The relational expression of each variable is defined below. Signal 4th line: Define the equation of dx / dt = -x

2. Compile the model into FMUs

Start JModelica in the model file directory. (Change the installation location as appropriate)

bash

/home/ubuntu/JModelica/bin/jm_ipython.sh

ipython

from pymodelica import compile_fmu
hello_fmu = compile_fmu("HelloWorld","./ode_test.mo")

3. Load FMUs

ipython

from pyfmi import load_fmu
hello_model = load_fmu(hello_fmu)

4. Calculate the loaded FMUs

Calculate for 1 second

ipython

res = hello_model.simulate(final_time=1)

5. View results

The result of the state variable x can be accessed with res ["x "]. Graph the above analytical solutions together.

ipython

import numpy as np
from matplotlib import pyplot as plt
t = np.linspace(0,1,101)
x = np.exp(-t) 
plt.plot(t, x, label="$x=e^(-t)$")
plt.plot(res["time"],res["x"],"--",label="JModelica")
plt.legend()
plt.show()

Solve ordinary differential equations with Assimulo

What is Assimulo

• Python module that handles the integrator of JModelica (can be used in any environment where JModelica can be used)
• Convenient for solving ordinary differential equations
• For Anaconda people, install with conda install -c https://conda.binstar.org/chria assimulo

ipython

from assimulo.solvers import CVode
from assimulo.problem import Explicit_Problem
#Define a function that represents a differential equation
def ode_func(t,x):
    dxdt = -x[0]
    return np.array([dxdt])
#Define and calculate a model that includes explicit problems and integrators
exp_mod = Explicit_Problem(ode_func, 1) #The initial value of x is 1
exp_sim = CVode(exp_mod)
t1, x1 = exp_sim.simulate(1)#Calculation time 1 second
#Result plot
plt.plot(t, x, label="$x=\exp(-t)$")#Analytical solution calculated earlier
plt.plot(res["time"],res["x"],'--',label="JModelica")#Numerical solution of JModelica
plt.plot(t1,x1,'-.',label="assimulo")#Numerical solution of Assimulo
plt.legend()
plt.show()

Summary

• SymPy derived the analytical solution of the initial value problem of the first-order linear ordinary differential equation.
• Confirmed the basic usage of JModelica
• The appearance of the calculation result (numerical solution) is better than that of the analytical solution.
• Convenient because the calculation result can be handled in Python
• Numerically solved ordinary differential equations with Assimulo
• Convenient if an ordinary differential equation is established
• I don't know how to read Assimulo