Whilst Dymola is good at solving nonlinear equations in models automatically, this post has a look at how to create a *function* that solves a nonlinear equation. A useful tool when required, as well as a way to get further insight into how nonlinear equations are solved.

## Why solve a nonlinear equation as a function?

There are some cases where a function has to be used, and using a model is not an option. Functions contain an *algorithm* section; top-down variable allocations are defined, with no way to define a nonlinear system as an *equation* section would.

For example, the **Modelica.Media **classes and the road models in the VeSyMA libraries make use of replaceable functions. In some cases, these functions may contain nonlinear equations. A method to solve these equations is described below.

## How are nonlinear equations solved in Dymola?

The Modelica equations are rearranged to represent the nonlinear differential equations in this form:

**fn**(**x**) =** 0**

Where the function **fn**() has nonlinear characteristics.

This system of equations is solved using a root finding method similar to the Newton-Raphson method. Iterative variables, **x**, are iteratively updated based on the gradient of **fn**(**x**), until **fn**(**x**) converges on **0 **(i.e. the residuals converge to zero).

The gradient of **fn**(**x**) is the Jacobian matrix:

Figure 1: Jacobian matrix

In the Newton-Raphson method the following update equation is used:

Figure 2. The update equation

This equation describes how the inverse of the Jacobian matrix is used to update the iterative variables (refer to http://fourier.eng.hmc.edu/e176/lectures/NM/node21.html for more detail), a similar method is used in Dymola. See the post about symbolic manipulation for further details.

## A simple nonlinear model example

This is a simple nonlinear model example. The closest point on a curve is found to a given point, which will be used throughout the rest of the post.

model LineContact

Real r[2] “Point”;

Real c[2](start={0,0}) “Closest point on valley function to Point”;

equation

r[1] = time;

r[2] = lineFunction(r[1]);

c[2] = valley(c[1]);

// Closest point on the valley function is perpendicular to (c-r)

(c-r)*{1,valleyGradient(c[1])} = 0;

end LineContact;

In the model the closest point, **c**, on the *valley* function, is found to point **r**. The *valley* function in this case is just a quadratic function, the *valleyGradient* function is the partial derivative of the *valley* function, with respect to the input. The *lineFunction* is any function that does not intersect the *valley* function.

An animation of this model is below:

### Obtaining the nonlinear equation and the Jacobian

When this model is translated, it generates the nonlinear equation as due to this equation:

**(c-r)*{1,valleyGradient(c[1])} = 0;**

To obtain the Jacobian and the nonlinear equation information about this model from Dymola do the following:

- Select
**Generate listing of translated Modelica code in dsmodel.mof**from the Simulation>Setup on the Translation tab (see Figure 3) - Set
**Advanced.OutputModelicaCodeWithJacobians = true**in the Dymola command line - Translate the model

Figure 3. Selecting **Generate listing of translated Modelica code dsmodel.mof**

Now that the model is translated, the dsmodel.mof file is generated which contains the nonlinear equation and the Jacobian as in Figure 4.

Figure 4. Nonlinear equation in dsmodel.mof

In Figure 4, the residual of the nonlinear equation is:

**(c-r)*{1,valleyGradient(c[1])} = 0;**

The Jacobian of the residual is also supplied as **J[:,:]**.

### Creating a function that can solve a nonlinear equation

The residual calculation and the Jacobian from the flat Modelica in Figure 4 can be used in an iterative loop to solve the nonlinear equation in a function as below:

// Iterative loop

while i < n and abs(residualValue) > threshold loop

i := i + 1;

JValue[1, 1] := 1.0 + (c[2] – r[2])*valleyGradient_der(c[1], 1.0) +

valleyGradient(c[1])*valley_der(c[1], 1.0);

JValueInv[1, 1] := 1/JValue[1, 1];

c[1] := c[1] – scalar(residualValue*JValueInv);

c[2] := valley(c[1]);

residualValue := (c – r)*{1,valleyGradient(c[1])};

end while;

The number of times the loop can repeat is limited *n* to avoid an infinite loop. To reduce computation the loop exits when the residual gets less than *threshold*. The code calculates the Jacobian of the residual and updates the iterative variable using the equation in Figure 4. Then the residual is recalculated and the loop is repeated until an exit condition is met.

This function was tested and it gave identical results to that of the nonlinear equation generated by the ‘*A simple nonlinear equation model*‘.

## Disadvantages of putting the nonlinear equation into a function

The nonlinear function above has some disadvantages to that of the automatically generated nonlinear equation generated in Dymola:

- Less robust nonlinear solver
- Logging of nonlinear solver for debugging purposes has to be manually created
- On every solve of the nonlinear equation the iterative variable is solved from an initial value. In the model the last value of the iteration variable is used as the initial value, this typically reduces computation significantly
- Requires significant amount of time to create

Most of disadvantages above can be worked around by moving the code to an external function and refining the algorithms further. This approach is used by the VeSyMA road models.

## Conclusion

This post describes how a function with nonlinear equations can be created. This solves the situation where a function is required that contains a nonlinear system. Possibly the Modelica standard should be extended to include a class like functions that support nonlinear equations.

The full code used in this example is available to download here

**Written by: Garron Fish – Chief Engineer**

**Please get in touch if you have any questions or have got a topic in mind that you would like us to write about. You can submit your questions / topics via: Tech Blog Questions / Topic Suggestion**