Perturb the solver's solution of a system's states to better satisfy time-invariant solution relationships - MATLAB mdlProjection

Required

No

Languages

C, C++

Syntax

#define MDL_PROJECTIONvoid mdlProjection(SimStruct *S)

Arguments

S: SimStruct representing an S-Function block.

Description

这种方法适用于我们e with S-functions that model dynamic systems whose states satisfy time-invariant relationships, such as those resulting from mass or energy conservation or other physical laws. The Simulink^®engine invokes this method at each time step after the model's solver has computed the S-function's states for that time step. Typically, slight errors in the numerical solution of the states cause the solutions to fail to satisfy solution invariants exactly. YourmdlProjectionmethod can compensate for the errors by perturbing the states so that they more closely approximate solution invariants at the current time step. As a result, the numerical solution adheres more closely to the ideal solution as the simulation progresses, producing a more accurate overall simulation of the system modeled by your S-function.

YourmdlProjectionmethod's perturbations of system states must fall within the solution error tolerances specified by the model in which the S-function is embedded. Otherwise, the perturbations may invalidate the solver's solution. It is up to yourmdlProjectionmethod to ensure that the perturbations meet the error tolerances specified by the model. SeePerturbing a System's States Using a Solution Invariantfor a simple method for perturbing a system's states. The following articles describe more sophisticated perturbation methods that yourmdlProjectionmethod can use.

C.W. Gear, “Maintaining Solution Invariants in the Numerical Solution of ODEs,”Journal on Scientific and Statistical Computing, Vol. 7, No. 3, July 1986.
L.F. Shampine, “Conservation Laws and the Numerical Solution of ODEs I,”Computers and Mathematics with Applications, Vol. 12B, 1986, pp. 1287–1296.
L.F. Shampine, “Conservation Laws and the Numerical Solution of ODEs II,”Computers and Mathematics with Applications, Vol. 38, 1999, pp. 61–72.

Example

Perturbing a System's States Using a Solution Invariant

Here is a simple, Taylor-series-based approach to perturbing a system's states. Suppose your S-function models a dynamic system having a solution invariant, $g (X, t)$ , i.e., $g$ is a continuous, differentiable function of the system states, $X$ , and time, $t$ , whose value is constant with time. Then

$X_{n} ≅ X_{n}^{*} + J_{n}^{T} {(J_{n} J_{n}^{T})}^{- 1} R_{n}$

where

$X_{n}$ is the system's ideal state vector at the solver's current time step
$X_{n}^{*}$ is the approximate state vector computed by the solver at the current time step
$J_{n}$ is the Jacobian of the invariant function evaluated at the point in state space specified by the approximate state vector at the current time step:

$J_{n} = \frac{\partial g}{\partial X} (X_{n}^{*}, t_{n})$
$t_{n}$ is the time at the current time step
$R_{n}$ is the residual (difference) between the invariant function evaluated at $X_{n}$ and $X_{n}^{*}$ at the current time step:

$R_{n} = g (X_{n}, t_{n}) - g (X_{n}^{*}, t_{n})$

Note

The value of $g (X_{n}, t_{n})$ is the same at each time step and is known by definition.

Given a continuous, differentiable invariant function for the system that your S-function models, this formula allows your S-function'smdlProjectionmethod to compute a perturbation

$J_{n}^{T} {(J_{n} J_{n}^{T})}^{- 1} R_{n}$

of the solver's numerical solution, $X_{n}^{*}$ , that more closely matches the ideal solution, $X_{n}$ , keeping the S-function's solution from drifting from the ideal solution as the simulation progresses.

MATLAB Example

This example illustrates how the perturbation method outlined in the previous section can keep a model's numerical solution from drifting from the ideal solution as a simulation progresses. Consider the following model,mdlProjectionEx1:

The PredPrey block references an S-function,predprey_noproj.m, that uses the Lotka-Volterra equations

$\begin{array}{l} \dot{x} = a x (1 - y) \\ \dot{y} = - c y (1 - x) \end{array}$

to model predator-prey population dynamics, where $x (t)$ is the population density of the predators and $y (t)$ is the population density of prey. The ideal solution to the predator-prey ODEs satisfies the time-invariant function

$x^{- c} e^{c x} y^{- a} e^{a y} = d$

where $a$ , $c$ , and $d$ are constants. The S-function assumesa = 1,c = 2, andd = 121.85.

The Invariant Residual block in this model computes the residual between the invariant function evaluated along the system's ideal trajectory through state space and its simulated trajectory:

$R_{n} = d - x_{n}^{- c} e^{c x_{n}} y_{n}^{- a} e^{a y_{n}}$

where $x_{n}$ and $y_{n}$ are the values computed by the model's solver for the predator and prey population densities, respectively, at the current time step. Ideally, the residual should be zero throughout simulation of the model, but simulating the model reveals that the residual actually strays considerably from zero:

Now consider the following model,mdlProjectionEx2:

This model is the same as the previous model, except that its S-function,predprey.m, includes amdlProjectionmethod that uses the perturbation approach outlined inPerturbing a System's States Using a Solution Invariant为了弥补数值漂移。因此,the numerical solution more closely tracks the ideal solution as the simulation progresses as demonstrated by the residual signal, which remains near or at zero throughout the simulation: