2. Initialization of JMUs

2.1. Solving DAE initialization problems

Before a model can be simulated it must be initialized, i.e. consistent initial values must be computed. To do this, JModelica.org supplies the JMUModel member function initialize, which initializes the JMUModel. The function is called after compiling and creating a JMUModel:

# Compile the stationary initialization model into a JMU
from pymodelica import compile_jmu
model_name = compile_jmu("My.Model", "/path/to/MyModel.mo")

# Load the model instance into Python
from pyjmi import JMUModel
init_model = JMUModel(model_name)
# Solve the DAE initialization system
init_result = init_model.initialize()

The JMUModel instance init_model is now initialized and is ready to be simulated.

The interactive help for the initialize method is shown by the command:

>>> help(init_model.initialize)   
    The initialization method depends on which algorithm is used, this can 
    be set with the function argument 'algorithm'. Options for the algorithm 
    are passed as option classes or as pure dicts. See 
    JMUModel.initialize_options for more details.
    The default algorithm for this function is IpoptInitializationAlg. 
        algorithm --
            The algorithm which will be used for the initialization is 
            specified by passing the algorithm class as string or class 
            object in this argument. 'algorithm' can be any class which 
            implements the abstract class AlgorithmBase (found in 
            algorithm_drivers.py). In this way it is possible to write own 
            algorithms and use them with this function.
            Default: 'IpoptInitializationAlg'
        options -- 
            The options that should be used in the algorithm. For details on 
            the options do:
                >> myModel = JMUModel(...)
                >> opts = myModel.initialize_options()
                >> opts?
            Valid values are: 
                - A dict which gives IpoptInitializationAlgOptions with 
                  default values on all options except the ones listed in 
                  the dict. Empty dict will thus give all options with 
                  default values.
                - An options object.
            Default: Empty dict
        Result object, subclass of algorithm_drivers.ResultBase.

Options for the available initialization algorithms can be set by first retrieving an options object using the JMUModel method initialize_options:

>>> help(init_model.initialize_options)
    Get an instance of the initialize options class, prefilled with default 
    values. If called without argument then the options class for the 
    default initialization algorithm will be returned.
        algorithm --
            The algorithm for which the options class should be fetched. 
            Possible values are: 'IpoptInitializationAlg', 'KInitSolveAlg'.
            Default: 'IpoptInitializationAlg'
        Options class for the algorithm specified with default values.

Having solved the initialization problem, the result of the initialization can be retrieved from the return result object:

x = init_result['x']
y = init_result['y']

2.2. How JModelica.org creates the initialization system of equations

To find a set of consistent initial values a system of non-linear equations, called the system of initialization equations, is solved. This system is composed from the DAE equations, the initial equations, some resulting from start attributes with the fixed attribute set to true. Start attributes with the fixed attribute set to false are treated as initial guesses for the numerical algorithm used to solve the initialization problem

Some initialization algorithms require the system of initial equations to be well defined in the sense that the number of variables must be equal to the number of equations. If this is not the case, the

  • If the number of equations is greater than the number of variables the system is overdetermined. Such a system may not have a solution, and will be treated as ill-defined. An exception is thrown in this case.

  • If the number of equations is less than the number of variables the system is underdetermined and such a system has infinitely many solutions. In this case, the compiler tries to balance the system by setting some fixed attributes to true. So if the user supplies too few initial conditions, some variables with the attribute fixed set to false may be changed to true during initialization.

2.3. Initialization algorithms

2.3.1. Initialization using IPOPT

JModelica.org provides a method for DAE initialization that is based on IPOPT, the mathematical formulation of the algorithm can be found in the JMI API documentation. Note that this algorithm does not rely on the causalization procedure (in particular the BLT transformation) which is common. Instead, the DAE residual is introduced as an equality constraint when solving an optimization problem where the squared difference between the non-fixed start values and their corresponding variables are minimized. As a consequence, the algorithm relies on decent start values for all variables. This approach is generally more sensitive to lacking initial guesses for start values than are algorithms based on causalization.

The algorithm provides the options summarized in Table 6.1.

Table 6.1. Options for the collocation-based optimization algorithm

Option Default Description
result_file_name Empty string (default generated file name will be used) Specifies the name of the file where the optimization result is written. Setting this option to an empty string results in a default file name that is based on the name of the optimization class.
result_format 'txt' Specifies in which format to write the result. Currently only textual mode is supported.
write_scaled_result False Write the scaled optimization result if set to True. This option is only applicable when automatic variable scaling is enabled. Only for debugging use.

In addition to the options for the collocation algorithm, IPOPT options can also be set by modifying the dictionary IPOPT_options contained in the collocation algorithm options object. Here, all valid IPOPT options can be specified, see the IPOPT documentation for further information. For example, setting the option max_iter:

opts['IPOPT_options']['max_iter'] = 300

makes IPOPT terminate after 300 iterations even if no optimal solution has been found.

Some statistics from IPOPT can be obtained by issuing the command:

>>> res_init.solver.init_opt_ipopt_get_statistics()

The return argument of this function can be found by using the interactive help:

>>> help(res_init.solver.init_opt_ipopt_get_statistics)
    Get statistics from the last optimization run.
        return_status -- 
            The return status from IPOPT.
        nbr_iter -- 
            The number of iterations. 
        objective -- 
            The final value of the objective function.
        total_exec_time -- 
            The execution time.

2.3.2. Initialization using KInitSolveAlg

JModelica.org also provides a method for DAE initialization based on the non-linear equation solver KINSOL from the SUNDIALS suite. KINSOL is currently comprised in the Assimulo package, included when installing JModelica.org. KINSOL is based on Newtons method for solving non-linear equations and is thus locally convergent. Attempts are made to make KInitSolveAlg as robust as possible but the possibility of finding a local minimum instead of the solution still remains. If the solution found by KInitSolveAlg is a local minimum a warning will be printed. The initial guesses passed to KINSOL are the ones supplied as start attributes in the current Modelica model.

KInitSolveAlg also implements an improved linear solver connected to KINSOL. This linear solver implements Tikhonov regularization to handle the problems of singular Jacobians as well as support for SuperLU, an efficient sparse linear solver.

The options providable are summarized in Table 6.2.

Table 6.2. Options for KInitSolveAlg

Option Default Description
use_constraints False A flag indicating whether constraints are to be used during initialization. Further explained in the section called “The use of constraints”.
constraints None A numpy.array containing floats that, when supplied, defines the constraints on the variables. Further explained in the section called “The use of constraints”.
result_format 'txt' Specifies in which format to write the result. Currently only textual mode is supported.
result_file_name Empty string (default generated file name will be used) Specifies the name of the file where the optimization result is written. Setting this option to an empty string results in a default file name that is based on the name of the optimization class.
KINSOL_options A dictionary with the defalt KINSOL options These are the options sent to the KINSOL solver. These are reviewed in detail in Table 6.3.

Table 6.3. Options for KINSOL contained in the KINSOL_options dictionary

Options Default Descriptions
use_jac True Flag indicating whether or not KINSOL uses the jacobian supplied by JModelica.org (True) or if KINSOL evaluates the Jacobian through finite differences (False). Finite differences is currently not available in sparse mode.
sparse False Flag indicating whether the problem should be treated as sparse (True) or dense (False).

The use of constraints

KINSOL, and hence also KInitSolvAlg, only supports simple unilateral constraints, that is constraining a variable to being positive or negative. If the option use_constraints is set to True, constraints are used. Which constraints that are used depends on whether or not the user has supplied constraints with the constraints option. If set, these will be used otherwise constraints will be computed by reading the min and max attributes from the Modelica file. How the constraint array is written is summarized in Table 6.4.

Table 6.4. Values allowed in the constraints array

Value Constraint
0.0 Unconstrained.
1.0 Greater than, or equal to, zero.
2.0 Greater than zero.
-1.0 Less than, or equal to, zero.
-2.0 Less than zero.

When the constraints are read from the Modelica file the value from Table 6.4 most fitting to the min and max values is chosen. For example a variable with min set to 3.2 and max set to 5.6 is constrained to be greater than zero. When the algorithm is finished however the result will be compared with the min and max values from the model testing if the solution fulfills the constraints set by the Modelica file.

Verbosity of KINSOL

There are four different levels of verbosity in KINSOL with 0 being silent and 3 being the most verbose. The verbosity level is controlled by the FMU log level. Table 6.5 describes what is output.

Table 6.5. Verbosity levels in KINSOL

FMU log level Verbosity level Output
<= 2 0 No information displayed.
3 1 In each nonlinear iteration the following information is displayed: the scaled Euclidean norm of the residual at the current iterate, the scaled Euclidian norm of the Newton step as well as the number of function evaluations performed so far.
4 2 Level 1 output as well as the Euclidian and in finity norm of the scaled residual at the current iterate
>= 5 3 Level 2 output plus additional values used by the global strategy as well as statistical information from the linear solver.