overall_system_ode_for_dha.m

一个matlab的将军模型

function xdot = overall_system_ode_for_DHA(t,x,flag,ode_param,p,x0)

% Get overall system derivative vector for the given the FSM state vector.
%
% Syntax:
%   "xdot = overall_system_ode(t,x,flag,ode_param)"
%
% Description:
%   Return the composite derivative vector "xdot" for the given FSM state
%   vector "q" and the SCSBs information are stored in "SCSBlocks".
%
% Implementation:
%   For the "k"-th SCSB, compute the vector "uk", the outputs of FSMBs in
%   "q" in the order that they feed into the "k"-th SCSB in the Simulink
%   diagram. Use "uk" to call the `switching function` for the SCSB to
%   obtain the derivative vector for that block.  Stack the derivative
%   vectors from the SCSBs together to form the overall (composite)
%   derivative vector.
%
% See Also:
%   piha,overall_system_clock,overall_system_matrix

global GLOBAL_PIHA
global evalnumber

evalnumber=evalnumber+1;
if nargin<=4
    p=[];
end;


q = ode_param;
xdot = [];
startx = 0;
startp = 0;
for k = 1:length(GLOBAL_PIHA.SCSBlocks)
  nx = GLOBAL_PIHA.SCSBlocks{k}.nx;
  nz = GLOBAL_PIHA.SCSBlocks{k}.nz;
  nup = GLOBAL_PIHA.SCSBlocks{k}.nup;
  np = GLOBAL_PIHA.SCSBlocks{k}.paradim;
  swfunck = GLOBAL_PIHA.SCSBlocks{k}.swfunc;
  uk = q(GLOBAL_PIHA.SCSBlocks{k}.fsmbindices);
  xdot = [xdot; deriv(swfunck,x(startx+1:startx + nx+nz),nx,nz,nup,uk,t,p(startp+1:startp + np),x0)];
  startx = startx + nx;
  startp = startp + np;
end


return

% -----------------------------------------------------------------------------

function xdot = deriv(swfunc,x,nx,nz,nup,u,t0,p,x0)

if (nargin(swfunc)>2 & ~isempty(p))
    
    %Calculate controller output for this period.  It should only be a function of the 
    %controller state (which does not change over a sample period) and the initial condition
    %of the plant state.
    [sys_temp,type] = feval(swfunc,[x0;zeros(nup,1)],u,p);

    %Time varying case
    %[sys_temp,type] = feval(swfunc,[x0;zeros(nup,1)],u,p,t0);
    
    [sys,type] = feval(swfunc,[x(1:nx+nz,1) ; sys_temp(nx+nz+1:nx+nz+nup,1) ],u,p);
    %Time varying case
    [sys,type] = feval(swfunc,[x(1:nx+nz,1) ; sys_temp(nx+nz+1:nx+nz+nup,1) ],u,p,t0);
    
    %[sys,type] = feval(swfunc,x,u,p);
else
    %Calculate controller output for this period.  It should only be a function of the 
    %controller state (which does not change over a sample period) and the initial condition
    %of the plant state.
    [sys_temp,type] = feval(swfunc,[x0;zeros(nup,1)],u);
    %Time varying case
    %[sys_temp,type] = feval(swfunc,[x0;zeros(nup,1)],u,t0);
   
    [sys,type] = feval(swfunc,[x(1:nx+nz,1) ; sys_temp(nx+nz+1:nx+nz+nup,1) ],u);
    %Time varying case
    %[sys,type] = feval(swfunc,[x(1:nx+nz,1) ; sys_temp(nx+nz+1:nx+nz+nup,1) ],u,t0);
    
    %[sys,type] = feval(swfunc,x,u);
end;

switch type
  case {'clock','nonlinear'}
    xdot(1:nx,1) = sys(1:nx,1);
    xdot(nx+1:nx+nz,1) = zeros(nz,1);
  case 'linear',
    % If the differential equation type is linear dx/dt = Ax + b, the returned
    % sys from swfunc is a structure containing the matrix A and the vector b.
    A = sys.A; b = sys.b;
    xdot = A*x + b;
  otherwise,
    error(['Unknown dynamics type ''' type '''!!!'])
end
return