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
