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

No

C, C++

```
#define MDL_PROJECTION
void mdlProjection(SimStruct *S)
```

`S`

SimStruct representing an S-Function block.

This method is intended for use 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. Your `mdlProjection`

method 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.

Your `mdlProjection`

method'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 your `mdlProjection`

method
to ensure that the perturbations meet the error tolerances specified
by the model. See Perturbing a System's States Using a Solution Invariant for
a simple method for perturbing a system's states. The following articles
describe more sophisticated perturbation methods that your `mdlProjection`

method
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.

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}\cong {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's `mdlProjection`

method
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.

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}=ax(1-y)\\ \dot{y}=-cy(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}^{cx}{y}^{-a}{e}^{ay}=d$$

where $$a$$, $$c$$, and $$d$$ are constants. The S-function
assumes `a = 1`

,` c = 2`

, and ```
d
= 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 a `mdlProjection`

method that uses the perturbation
approach outlined in Perturbing a System's States Using a Solution Invariant
to compensate for numerical drift. As a result, 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: