During the initialization phase of a dynamic simulation problem, it often happens that large nonlinear systems of equations must be solved by means of an iterative solver. The convergence of such solvers critically depends on the choice of initial guesses for the unknown variables. The process can be made more robust by providing an alternative, simplified version of the model, such that convergence is possible even without accurate initial guess values, and then by continuously transforming the simplified model into the actual model. This transformation can be formulated using expressions of this kind:

lambda*actual + (1-lambda)*simplified

in the formulation of the system equations, and is usually called a homotopy transformation. If the simplified expression is chosen carefully, the solution of the problem changes continuously with lambda, so by taking small enough steps it is possible to eventually obtain the solution of the actual problem.

It is recommended to perform (conceptually) one homotopy iteration over the whole model, and not several homotopy iterations over the respective non-linear algebraic equation systems. The reason is that the following structure can be present:

w = f// has homotopy operator_{1}(x)0 = f_{2}(der(x), x, z, w)

Here, a non-linear equation system
** f** is present. The
homotopy operator is, however used on a variable that is an 'input'
to the non-linear algebraic equation system, and modifies the
characteristics of the non-linear algebraic equation system. The
only useful way is to perform the homotopy iteration over

`f`_{1}

`f`_{2}

The suggested approach is 'conceptual', because more efficient implementations are possible, e.g., by determining the smallest iteration loop, that contains the equations of the first BLT block in which a homotopy operator is present and all equations up to the last BLT block that describes a non-linear algebraic equation system.

A trivial implementation of the homotopy operator is obtained by defining the following function in the globalscope:

functionhomotopyinputReal actual;inputReal simplified;outputReal y;algorithmy := actual;annotation(Inline = true);endhomotopy;

homotopy(actual=actual, simplified=simplified)

The scalar expressions 'actual' and 'simplified' are subtypes of Real. A Modelica translator should map this operator into either of the two forms:

- Returns 'actual'
*[a trivial implementation]*. -
In order to solve algebraic systems of equations, the operator might during the solution process return a combination of the two arguments, ending at actual,

*e.g.,**actual*lambda + simplified*(1-lambda),*where

`lambda`

is a homotopy parameter going from 0 to 1.The solution must fulfill the equations for homotopy returning 'actual'.

In electrical systems it is often difficult to solve non-linear algebraic equations if switches are part of the algebraic loop. An idealized diode model might be implemented in the following way, by starting with a 'flat' diode characteristic and then move with the homotopy operator to the desired 'steep' characteristic:

modelIdealDiode ...parameterReal Goff = 1e-5;protectedRealGoff_flat = max(0.01, Goff);RealGoff2;equationoff = s < 0; Goff2 =homotopy(actual=Goff, simplified=Goff_flat); u = s*(ifoffthen1elseRon2) + Vknee; i = s*(ifoffthenGoff2else1 ) + Goff2*Vknee; ...endIdealDiode;

In electrical systems it is often useful that all voltage sources start with zero voltage and all current sources with zero current, since steady state initialization with zero sources can be easily obtained. A typical voltage source would then be defined as:

modelConstantVoltageSourceextendsModelica.Electrical.Analog.Interfaces.OnePort;parameterModelica.SIunits.Voltage V;equationv =homotopy(actual=V, simplified=0.0);endConstantVoltageSource;

In fluid system modelling, the pressure/flowrate relationships are highly nonlinear due to the quadratic terms and due to the dependency on fluid properties. A simplified linear model, tuned on the nominal operating point, can be used to make the overall model less nonlinear and thus easier to solve without accurate start values. Named arguments are used here in order to further improve the readability.

modelPressureLossimportSI = Modelica.SIunits; ...parameterSI.MassFlowRate m_flow_nominal "Nominal mass flow rate";parameterSI.Pressure dp_nominal "Nominal pressure drop"; SI.Density rho "Upstream density"; SI.DynamicViscosity lambda "Upstream viscosity";equation... m_flow =homotopy(actual = turbulentFlow_dp(dp, rho, lambda), simplified = dp/dp_nominal*m_flow_nominal); ...endPressureLoss;

Note that the homotopy operator **shall not** be
used to combine unrelated expressions, since this can generate
singular systems from combining two well-defined systems.

modelDoNotUse Real x;parameterReal x0 = 0;equationder(x) = 1-x;initial equation0 =homotopy(der(x), x - x0);endDoNotUse;

The initial equation is expanded into

0 = lambda*der(x) + (1-lambda)*(x-x0)

and you can solve the two equations to give

x = (lambda+(lambda-1)*x0)/(2*lambda - 1)

which has the correct value of `x0`

at ```
lambda =
0
```

and of `1`

at `lambda = 1`

, but
unfortunately has a singularity at `lambda = 0.5`

.

Generated at 2019-07-19T01:39:52Z by OpenModelicaOpenModelica 1.14.0~dev-26610-gea7fc62 using GenerateDoc.mos