.ModelicaReference.Operators.'homotopy()'

Information

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 = f1(x) // has homotopy operator
0 = f2(der(x), x, z, w)

Here, a non-linear equation system f2 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 f1 and f2 together.

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:

function homotopy
  input Real actual;
  input Real simplified;
  output Real y;
algorithm
  y := actual;
  annotation(Inline = true);
end homotopy;

Syntax

homotopy(actual=actual, simplified=simplified)

Description

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

1. Returns 'actual' [a trivial implementation].

2. 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'.

Examples

Example 1

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:

model IdealDiode
  ...
  parameter Real Goff = 1e-5;
protected
  Real Goff_flat = max(0.01, Goff);
  Real Goff2;
equation
  off = s < 0;
  Goff2 = homotopy(actual=Goff, simplified=Goff_flat);
  u = s*(if off then 1     else Ron2) + Vknee;
  i = s*(if off then Goff2 else 1   ) + Goff2*Vknee;
  ...
end IdealDiode;

Example 2

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:

model ConstantVoltageSource
  extends Modelica.Electrical.Analog.Interfaces.OnePort;
  parameter Modelica.Units.SI.Voltage V;
equation
  v = homotopy(actual=V, simplified=0.0);
end ConstantVoltageSource;

Example 3

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.

model PressureLoss
  import SI = Modelica.Units.SI;
  ...
  parameter SI.MassFlowRate m_flow_nominal "Nominal mass flow rate";
  parameter SI.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);
  ...
end PressureLoss;

Example 4

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.

model DoNotUse
  Real x;
  parameter Real x0 = 0;
equation
  der(x) = 1-x;
initial equation
  0 = homotopy(der(x), x - x0);
end DoNotUse;

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.