This blog post looks at how partial derivatives of a model can be calculated by exporting a model as an FMU, then importing the FMU into Dymola and using the fmiGetDirectionalDerivative() function.
How can partial derivatives be calculated?
Modelica does not supply an easy way to calculate partial derivatives of a model or system. There are ways for calculating partial derivatives of functions, described in Using Automatic Differentiation for Partial Derivatives of Functions in Modelica, and the PDE library, which calculates partial derivatives for the use of creating models that use partial derivatives.
Fortunately, the FMI standard defines a fmi2GetDirectionalDerivative() function for calculating partial derivatives (see page 26 of the FMI-Specification). Therefore, one method to get the partial derivative of a model is to create an FMU of the model, and then import the FMU into Dymola and use this partial derivative function.
Calculating partial derivatives using an FMU
Creating and importing an FMU
For example, let’s take a simple driven pendulum model, Figure 1, and export it.
Figure 1. A simple driven pendulum model
The driven pendulum model can be exported as an FMU by pressing Ctrl+F9; see this blog post for more details. It can be imported by going to File>Open>Import FMU, either as “Model Exchange” or “Co-simulate” as both support the calculation of partial derivatives.
Modifying the FMU so interfacing functions can be used
When Dymola imports FMUs, it creates all the functions for interfacing with the FMU dll within the FMU code in Dymola. These FMU functions include the fmiGetDirectionalDerivative() function for calculating partial derivatives. However, these functions are stored in a protected section. The FMU model has to be modified to be able to make use of these functions as in Figure 2.
Figure 2. Commenting out protected so FMU interfacing functions can be used
Comment out the protected statement in the Dymola FMU code (as in Figure 2) that is found above the following line:
constant Integer fmi_NumberOfEventIndicators = 0;
This will allow the FMU interfacing functions to be used.
The fmiGetDirectionalDerivative() function
The fmiGetDirectionalDerivative() function in the Dymola FMU is used in the same way the fmi2GetDirectionalDerivative() function described on page 26 of the FMI-Specification. Except the nUnknown and nKnown inputs are not present in the Dymola fmiGetDirectionalDerivative() function, as in Figure 3.
Figure 3. Inputs and outputs of the fmiGetDirectionalDerivative() function
The fmi input is the index to the external fmi object (see the external object blog). The z_refs and v_refs inputs, in Figure 3, are the valueReference values of the signals. The valueReference for rev.w is obtained from the modelDescription.xml file below:
<ScalarVariable
name="rev.w"
valueReference="33554433"
description="First derivative of angle phi (relative angular velocity)"
initial="exact">
<Real
declaredType="Modelica.Units.SI.AngularVelocity"
start="0"
derivative="356"/>
</ScalarVariable>The valueReference values are the indexes to the signals in the FMU.
The z_ref signal is the output or derivative signal for which the partial derivative is to be calculated, taken with respect to the v_ref signal which is the input or state signal. dv is a gain applied to the partial derivative. Please refer to page 26 of the FMI-Specification for a complete description.
Using the fmiGetDirectionalDerivative() function
I recommend creating a separate model, placing the modified FMU (see Figure 2) and adding in a Modelica.Blocks.Source.RealExpression block to this model as in Figure 4.
Figure 4. FMU block with partial derivatives used to calculate A and B matrices
Then set the real expression component as:
<fmu component name>.fmi_Functions.fmiGetDirectionalDerivative(<fmu component name>.fmi, z_ref, v_ref, {1}))Where z_ref and v_ref are the valueReference signal vectors as described in “The fmiGetDirectionalDerivative() function” section, above. In Figure 4 the Real expressions components are matrices and the fmiGetDirectionalDerivative() is called multiple times to populate these matrices.
How can partial derivatives be used?
Partial derivatives of models can have many uses, two cases where they are used is for calculating state observers such as in the extended Kalman filter and for control purposes such as LQR.
A nonlinear LQR controller example
A LQR controller is designed to be used with linear systems. This example looks at a way of using LQR controllers with nonlinear models by calculating the A and B matrices at each instance using fmiGetDirectionalDerivative() and then updating the LQR controller accordingly. This example is purely for demonstration purposes.
FMU plus partial derivatives model
A model of the driven pendulum in Figure 1 was exported as an FMU and imported into Dymola. This FMU was modified as in Figure 2 to allow access to fmiGetDirectionalDerivative(). A model with partial derivatives was created as in Figure 4.
Calculate LQR K value
The infinite-horizon continuous time LQR method was implemented as described in the wikipedia entry. This required solving P in the Riccati equation:
Figure 5. The Ricatti equation. N was set to zero in this example. Figure: wikipedia.
Solving the above equation introduced four nonlinear equations for this example. In some cases better results were obtained by supplying an initial guess value for P.
Simulation of the driven pendulum
The full controlled driven pendulum model was created as in Figure 6.
Figure 6. Full driven pendulum model with lqr controller that is updated
The FMI value calculates the A and B matrices which are passed to the LQR designer which determines the value of K. This K value is negated and multiplied with the states of the system to determine controlling torque of the driven pendulum in Figure 6.
The revolute angle is 0 (radians) when the pendulum is inverted and is pi at the bottom. In this example the pendulum is initially stationary at the bottom and then is driven to the inverted position.
The value for Q (see LQR wiki) had to be set to Q=[100,0;0,1], low values for Q[1,1] would result in the pendulum getting “stuck” while going up. The R value was set to 1.
Figure 7. Driven pendulum simulated using nonlinear LQR
Conclusion
This post describes how to calculate partial derivatives of a Modelica model and an example is provided on how these partial derivatives can be used to control a driven pendulum.
The code for the driven pendulum is available on request, via: Driven Pendulum Code
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
