Notebook

モデルを作ってシミュレーションを実行する¶

JModelica Users Guid - Version 2.10 の 5.2. A first example に簡単なモデルを作ってシミュレーションを実行し、結果をプロットする方法が書かれています。これを実行してみます。

Van der Pol oscillator の支配方程式を次に示します。

$\frac{d^2 x}{dt^2} - \mu (1-x^2) \frac{d x}{d t} + x = 0$

$x=x_2, \frac{dx}{dt} = x_1, \mu = 1$ とおくと次のように変形できます。

$\frac{d x_1}{d t} = (1-x_2^2) x_1 - x_2, \ \frac{d x_2}{dt} = x_1$ これは、ここで解く方程式

$\frac{d x_1}{d t} = (1-x_2^2) x_1 - x_2 + u, \ \frac{d x_2}{dt} = x_1$ の制御信号(control signal) を

$u=0$ としたものになります。初期条件を

$x_1(0)=0, \ x_2(0)=1$ とします。

モデルのソースコードを作成します。短いソースコードのモデルは Jupyter Notebook の Magicコマンド %%writefile で作成できます。

In [1]:

%%writefile VDP.mo
model VDP
    // State start values
    parameter Real x1_0 = 0;
    parameter Real x2_0 = 1;
    
    // The states
    Real x1(start = x1_0);
    Real x2(start = x2_0);
    
    // The control signal
    input Real u;
    
equation
    der(x1) = (1 - x2^2) * x1 - x2 + u;
    der(x2) = x1;
end VDP;

Overwriting VDP.mo

モデルをコンパイルして FMU を作成してロードします。

In [2]:

# Import the function for compilation of models and the load_fmu method
from pymodelica import compile_fmu
from pyfmi import load_fmu
# Compile model
fmu_name = compile_fmu("VDP","VDP.mo")
# Load model
vdp = load_fmu(fmu_name)

シミュレーションを実行します。

In [3]:

res = vdp.simulate(final_time=10)

Final Run Statistics: --- 

 Number of steps                                 : 148
 Number of function evaluations                  : 208
 Number of Jacobian evaluations                  : 3
 Number of function eval. due to Jacobian eval.  : 6
 Number of error test failures                   : 11
 Number of nonlinear iterations                  : 204
 Number of nonlinear convergence failures        : 0

Solver options:

 Solver                   : CVode
 Linear multistep method  : BDF
 Nonlinear solver         : Newton
 Linear solver type       : DENSE
 Maximal order            : 5
 Tolerances (absolute)    : 1e-06
 Tolerances (relative)    : 0.0001

Simulation interval    : 0.0 - 10.0 seconds.
Elapsed simulation time: 0.0104110240936 seconds.

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

In [4]:

%matplotlib notebook
# Import the plotting library
import matplotlib.pyplot as plt 
x1 = res['x1']
x2 = res['x2']
t = res['time']
plt.figure(1)
plt.plot(t, x1, t, x2)
plt.legend(('x1','x2'))
plt.title('Van der Pol oscillator.')
plt.ylabel('Angle (rad)')
plt.xlabel('Time (s)')
plt.show()

シミュレーション結果に含まれる変数や定数、パラメータは以下のようにして調べることができます。

In [5]:

for name in res.keys():
    print(name)

time
x1_0
x2_0
_block_jacobian_check_tol
_cs_rel_tol
_cs_step_size
_events_default_tol
_events_tol_factor
_nle_jacobian_finite_difference_delta
_nle_solver_default_tol
_nle_solver_max_residual_scaling_factor
_nle_solver_min_residual_scaling_factor
_nle_solver_min_tol
_nle_solver_regularization_tolerance
_nle_solver_step_limit_factor
_nle_solver_tol_factor
_time_events_default_tol
der(x1)
der(x2)
x1
x2
u
_block_solver_experimental_mode
_cs_experimental_mode
_cs_solver
_iteration_variable_scaling
_log_level
_nle_active_bounds_mode
_nle_jacobian_calculation_mode
_nle_jacobian_update_mode
_nle_solver_exit_criterion
_nle_solver_max_iter
_nle_solver_max_iter_no_jacobian
_residual_equation_scaling
_block_jacobian_check
_block_solver_profiling
_enforce_bounds
_nle_brent_ignore_error
_nle_solver_check_jac_cond
_nle_solver_use_last_integrator_step
_nle_solver_use_nominals_as_fallback
_rescale_after_singular_jac
_rescale_each_step
_runtime_log_to_file
_use_Brent_in_1d
_use_jacobian_equilibration
_use_newton_for_brent

最初がアンダーバーで始まるものを除くと以下のようになります。

In [6]:

for name in res.keys():
    if name[0] != '_':
        print (name)

time
x1_0
x2_0
der(x1)
der(x2)
x1
x2
u