mdlProjection

Возмущение решения решателем состояний системы для лучшего удовлетворения отношений инвариантного по времени решения

Необходимый

Нет

Языки

C, C++

Синтаксис

#define MDL_PROJECTION
void mdlProjection(SimStruct *S)

Аргументы

S: SimStruct представляет S-функциональный блок.

Описание

Этот способ предназначен для использования с S-функциями, которые моделируют динамические системы, состояния которых удовлетворяют временным инвариантным соотношениям, таким как те, которые являются результатом сохранения массы или энергии или других физических законов. Модуль Simulink ® вызывает этот метод на каждом временном шаге после того, как решатель модели вычислит состояния S-функции для этого временного шага. Обычно небольшие ошибки в числовом решении состояний приводят к тому, что решения не удовлетворяют инвариантам решения в точности. ВашmdlProjection способ может компенсировать ошибки путем возмущения состояний так, чтобы они более близко приближали инварианты решения на текущем временном шаге. В результате численное решение по мере продвижения моделирования более тесно согласуется с идеальным решением, обеспечивая более точное общее моделирование системы, смоделированной с помощью S-функции.

Ваш mdlProjection возмущения системного состояния метода должны находиться в пределах допусков ошибок решения, заданных моделью, в которую встроена S-функция. В противном случае возмущения могут сделать решение решателя недействительным. Это зависит от вашего mdlProjection для обеспечения соответствия возмущений допускам ошибок, заданным моделью. Простой метод возмущения состояний системы см. в разделе Возмущение состояний системы с помощью инварианта решения. В следующих статьях описаны более сложные методы возмущения, которые mdlProjection способ может быть использован.

C.W. Gear, «Поддержание инвариантов решения в числовом решении ОДУ», Journal on Scientific and Statistical Computing, Vol. 7, No. 3, July 1986.
Л. Ф. Шампин, «Законы сохранения и числовое решение ОДУ I», Компьютеры и математика с приложениями, т. 12B, 1986, стр. 1287-1296.
Л. Ф. Шампин, «Законы сохранения и числовое решение ОДУ II», Computers and Mathematics with Applications, Vol. 38, 1999, pp. 61-72.

Пример

Нарушение состояния системы с помощью инварианта решения

Вот простой, основанный на серии Тейлора подход к возмущению состояний системы. Предположим, что ваша S-функция моделирует динамическую систему, имеющую инвариант решения, $g (X,$ t), т.е. g является непрерывной, дифференцируемой функцией состояний системы, X и времени $,$ t, значение которой является постоянным со временем. Тогда

$_{}_{}^{}_{}^{Xn≅Xn*+JnT} {(_{}_{}^{} JnJnT}^{)} -_{}$ 1Rn

где

$_{Xn}$ - вектор идеального состояния системы на текущем шаге времени решателя;
$_{Xn}^{}$ * - приближенный вектор состояния, вычисленный решателем на текущем шаге времени
$_{Jn}$ - якобиан инвариантной функции, оцениваемый в точке пространства состояний, заданной вектором приближенного состояния на текущем шаге времени:
$_{} \frac{}{Jn=∂g∂X} (_{}^{Xn}_{*,}$ tn)
$_{tn}$ - время на текущем шаге времени
$_{Rn}$ представляет собой остаток (разность) между инвариантной функцией, оцениваемой при $_{Xn}$ и $_{Xn}^{}$ * на текущем временном шаге:
$_{Rn} = g_{} (_{Xn}, tn)_{}^{-} g_{} ($ Xn *, tn)
Примечание
Величина $g (_{} Xn,_{}$ tn) одинакова на каждом временном этапе и известна по определению.

Учитывая непрерывную, дифференцируемую инвариантную функцию для системы, которую моделирует ваша S-функция, эта формула позволяет вашей S-функции mdlProjection способ вычисления возмущения

$_{}^{JnT} {(_{}_{}^{} JnJnT}^{)} -_{}$ 1Rn

числового решения решателя Xn $_{}^{}$ *, которое более точно соответствует идеальному решению Xn, $_{}$ удерживая решение S-функции от дрейфа от идеального решения по мере продвижения моделирования.

Пример MATLAB

Этот пример иллюстрирует, как метод возмущений, описанный в предыдущем разделе, может удерживать численное решение модели от дрейфа от идеального решения в процессе моделирования. Рассмотрим следующие model,mdlProjectionEx1:

Блок PredPrey ссылается на S-функцию, которая использует уравнения Лотка-Вольтерра:

predprey_noproj.m

function predprey_noproj(block)
% This is one of a pair S-functions that model a population of predators
% and preys, the other being predprey.m.

% Both S-functions use the Lotka-Volterra equations
%
%   xdot = ax(1-y)
%   ydot = -cy(1-x)
%
% to model changes in the population where x(t) and y(t) are 
% the number of predators and prey, respectively, and 
% a and c are constants. Both S-functions assume a=1 and c=2 and
% that the population consists initially of one predator and three prey.
% 
% The solution to the Lotka-Volterra equations obeys the time-invariant
% relationship
%
%   (x^-c)*exp(cx)*(y^-a)*exp(ay) = d
%
% where d is a constant.
%
% The other S-function uses this relationship to compensate for
% numerical errors in the solver's solution of the model's state
% equations. This S-function does not attempt to compensate for the
% numerical errors. As a result, the solution drifts from the ideal
% solution as the simulation progresses.
%
% For more information, see the documentation for the mdlProjection 
% method in the online Simulink documentation.
%   
% Copyright 2006 The MathWorks, Inc.

  setup(block);
  
%endfunction

function setup(block)
  
  %% Register number of input and output ports
  block.NumInputPorts  = 0;
  block.NumOutputPorts = 2;

  block.OutputPort(1).Dimensions       = 1;
  block.OutputPort(1).SamplingMode = 'Sample';
  
  block.OutputPort(2).Dimensions       = 1;
  block.OutputPort(2).SamplingMode = 'Sample';

  
  %% Set block sample time to variable sample time
  block.SampleTimes = [0 0];
  
  %% Setup Dwork
  block.NumContStates = 2;
  
  %% Register methods
  block.RegBlockMethod('InitializeConditions',    @InitConditions); 
  block.RegBlockMethod('Outputs',                 @Output);  
  block.RegBlockMethod('Derivatives', @Derivatives);
  block.RegBlockMethod('Projection', @Projection);
  block.RegBlockMethod('Update', @Update);
  
%endfunction

function InitConditions(block)

  %% Initialize Dwork: 1 predator, 3 prey.
  block.ContStates.Data = [1 3]';
  
  
%endfunction

function Output(block)

  % Output number of predators.
  block.OutputPort(1).Data = block.ContStates.Data(1);
  
  % Output number of prey.
  block.OutputPort(2).Data = block.ContStates.Data(2);
  
  block.ContStates.Data;
  
%endfunction

function Derivatives(block)

states = block.ContStates.Data;
% System ODEs
block.Derivatives.Data = [2*states(1)*(1-states(2)); -states(2)*(1-states(1))];
  
%endfunction

function Projection(block)
    %not implemented
%endfunction

function Update(block)

%endfunction

$\begin{array}{l} \overset{}{} x˙=ax (1 - \\ \overset{y}{}) y˙=−cy \end{array}$ (1 − x)

для моделирования динамики популяции хищника-добычи, где $x (t$ ) - плотность популяции хищников, а $y ($ t) - плотность популяции добычи. Идеальное решение ODE хищника-добычи удовлетворяет инвариантной во времени функции

$x^{-}^{}^{cecxy}^{-}$ aeay = d

где $a$ , $c$ и $d$ - константы. S-функция предполагает a = 1, c = 2, и d = 121.85.

Блок Инвариантный остаток (Invariant Resident) в этой модели вычисляет остаток между инвариантной функцией, оцениваемой вдоль идеальной траектории системы через пространство состояний, и ее моделируемой траекторией:

$_{Rn} = d_{}^{-}^{xn_{-}}_{}^{}^{_{cecxnyn}}$ − aeayn

где $_{}$ xnand $_{}$ ynare значения, вычисленные решателем модели для плотности популяции хищника и добычи, соответственно, на текущем шаге времени. В идеале остаток должен быть равен нулю на протяжении всего моделирования модели, но моделирование модели показывает, что остаток на самом деле значительно отклоняется от нуля:

Рассмотрим следующую модель, mdlProjectionEx2:

Эта модель аналогична предыдущей, за исключением ее S-функции, predprey.m, включает mdlProjection метод, который использует подход возмущений, описанный в разделе Возмущение состояний системы с использованием инварианта решения для компенсации численного дрейфа. В результате численное решение более точно отслеживает идеальное решение по мере продвижения моделирования, как показано остаточным сигналом, который остается близким или равным нулю на протяжении всего моделирования:

predprey.m

function predprey(block)
% This is one of a pair S-functions that model a population of predators
% and preys, the other being predprey_noproj.m.

% Both S-functions use the Lotka-Volterra equations
%
%   xdot = ax(1-y)
%   ydot = -cy(1-x)
%
% to model changes in the population where x(t) and y(t) are 
% the number of predators and prey, respectively, and 
% a and c are constants. Both S-functions assume a=1 and c=2 and
% that the population consists initially of one predator and three prey.
% 
% The solution to the Lotka-Volterra equations obeys the time-invariant
% relationship
%
%   (x^-c)*exp(cx)*(y^-a)*exp(ay) = d
%
% where d is a constant.
%
% This S-function uses this relationship to compensate for
% numerical errors in the solver's solution of the model's state
% equations. The other S-function does not. As a result, 
% the other's solution drifts from the ideal solution as the 
% simulation progresses.
%
% For more information, see the documentation for the mdlProjection 
% method in the online Simulink documentation.
%   
% Copyright 2006 The MathWorks, Inc.
  setup(block);
  
%endfunction

function setup(block)
  
  %% Register number of input and output ports
  block.NumInputPorts  = 0;
  block.NumOutputPorts = 2;

  block.OutputPort(1).Dimensions       = 1;
  block.OutputPort(1).SamplingMode = 'Sample';
  
  block.OutputPort(2).Dimensions       = 1;
  block.OutputPort(2).SamplingMode = 'Sample';

  
  %% Set block sample time to variable sample time
  block.SampleTimes = [0 0];
  
  %% Setup Dwork
  block.NumContStates = 2;
  
  %% Register methods
  block.RegBlockMethod('InitializeConditions',    @InitConditions); 
  block.RegBlockMethod('Outputs',                 @Output);  
  block.RegBlockMethod('Derivatives', @Derivatives);
  block.RegBlockMethod('Projection', @Projection);
  
%endfunction

function InitConditions(block)

  %% Initialize Dwork: 1 predator, 3 prey.
  block.ContStates.Data = [1 3]';
  
  
%endfunction

function Output(block)

  % Output number of predators.
  block.OutputPort(1).Data = block.ContStates.Data(1);
  
  % Output number of prey.
  block.OutputPort(2).Data = block.ContStates.Data(2);
  block.ContStates.Data;
  
%endfunction

function Derivatives(block)

states = block.ContStates.Data;
% System ODEs
block.Derivatives.Data = [2*states(1)*(1-states(2)); -states(2)*(1-states(1))];
  
%endfunction

function Projection(block)

states = block.ContStates.Data;

%Computing the Jacobian of the Invariant
J = localInvJac(states);

%Computing the residual of the invariant
R = localInvRes(states);

%Perturbing the states
states = states + (J')*((J*(J'))^-1) * R;

%Update the states
block.ContStates.Data = states; 

%endfunction

function J = localInvJac(states)

%Computing the Jacobian of the invariant
L = 1/states(1) * exp(states(1));
R = 1/(states(2)^2)*exp(2*states(2));
Lx = (states(1)-1)/(states(1)^2)*exp(states(1));
Ry = 2*(states(2)-1)/(states(2)^3)*exp(2*states(2));

J = [Lx*R L*Ry];

%endfunction

function R = localInvRes(states)

%Invariant value
inv0 = 1.218481287142732e+002;

%Computing invariant
L = 1/states(1) * exp(states(1));
R = 1/(states(2)^2)*exp(2*states(2));
inv = L*R;

%Residual
R = inv0-inv;

%endfunction

См. также

Projection

Представлен в R2006b

Документация