clk_flow.m

一个matlab的将军模型

function R = clk_flow(X0,v,T)

% Transform a given set according to dynamic v and the time stamp T.
% 
% Syntax:
%   "R = clk_rch(X0,v,t)"
%
% Description:
%   Given the initial polytope "X0", compute the transformation of the set from "X0" 
%   in the clock direction given by the
%   clock vector "v" and time stamp T. 
%
% Implementation:


vertset = vertices(X0);
vertout = vertices;
for i=1:length(vertset)
   vertout
if isempty(dE)
  R = reach_no_eq(CI,dI,v);
else
  R = reach_one_eq(CE,dE,CI,dI,v);
end
return

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

function R = reach_no_eq(CI,dI,v)

% [CEreach,dEreach,CIreach,dIreach] are the constraints for the reachable 
% set.
CEreach = []; dEreach = [];
CIreach = []; dIreach = [];

% Go through the inequality contraints
Tplus = [];
Tminus = 0;
for i = 1:length(dI)
  if (CI(i,:)*v == 0)  % case (1.1)
    CIreach = [CIreach; CI(i,:)];
    dIreach = [dIreach; dI(i)];
  else
    if (CI(i,:)*v > 0) % case (1.2)
      Tminus = [Tminus i];
    else
      Tplus = [Tplus i];
    end
  end
end

% Case (1.2) (continue): compute all possible conjunctions of two
% inequalities in T+ and T-
for i = 1:length(Tminus)
  for j = 1:length(Tplus)
    if (Tminus(i) == 0)
      CIj = CI(Tplus(j),:);  dIj = dI(Tplus(j));
      c = -CIj/(CIj*v);
      d = -dIj/(CIj*v);
    else
      CIi = CI(Tminus(i),:); dIi = dI(Tminus(i));
      CIj = CI(Tplus(j),:);  dIj = dI(Tplus(j));
      c = CIi/(CIi*v) - CIj/(CIj*v);
      d = dIi/(CIi*v) - dIj/(CIj*v);
    end
    % normalize the normal vector
    magnitude = sqrt(c*c');
    if (magnitude == 0) 
      if (d < 0)
        error('Infeasible constraint found! Something is wrong!')
      end
    else
      c = c/magnitude;
      d = d/magnitude;
      CIreach = [CIreach; c];
      dIreach = [dIreach; d];
    end
  end
end

% create a linearcon object and return the reachable set
R = linearcon(CEreach,dEreach,CIreach,dIreach);
return

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

function R = reach_one_eq(CE,dE,CI,dI,v)

% [CEreach,dEreach,CIreach,dIreach] are the constraints for the reachable 
% set.
CEreach = []; dEreach = [];
CIreach = []; dIreach = [];

% Check the equality constraint
if (CE*v == 0)  % case (2.1.1)
  CEreach = [CEreach; CE];
  dEreach = [dEreach; dE];
  % Go through the inequality contraints
  Tplus = [];
  Tminus = 0;
  for i = 1:length(dI)
    if (CI(i,:)*v == 0)  % case (1.1)
      CIreach = [CIreach; CI(i,:)];
      dIreach = [dIreach; dI(i)];
    else
      if (CI(i,:)*v > 0) % case (1.2)
        Tminus = [Tminus i];
      else
        Tplus = [Tplus i];
      end
    end
  end
  % Case (1.2) (continue): compute all possible conjunctions of two
  % inequalities in T+ and T-
  for i = 1:length(Tminus)
    for j = 1:length(Tplus)
      if (Tminus(i) == 0)
        CIj = CI(Tplus(j),:);  dIj = dI(Tplus(j));
        c = -CIj/(CIj*v);
        d = -dIj/(CIj*v);
      else
        CIi = CI(Tminus(i),:); dIi = dI(Tminus(i));
        CIj = CI(Tplus(j),:);  dIj = dI(Tplus(j));
        c = CIi/(CIi*v) - CIj/(CIj*v);
        d = dIi/(CIi*v) - dIj/(CIj*v);
      end
      % normalize the normal vector
      magnitude = sqrt(c*c');
      if (magnitude == 0) 
        if (d < 0)
          error('Infeasible constraint found! Something is wrong!')
        end
      else
        c = c/magnitude;
        d = d/magnitude;
        CIreach = [CIreach; c];
        dIreach = [dIreach; d];
      end
    end
  end
else
  if (CE*v > 0) % case (2.1.2)
    CIreach = [CIreach; -CE];
    dIreach = [dIreach; -dE];
  else
    CIreach = [CIreach; CE];
    dIreach = [dIreach; dE];
  end
  % Go through inequality constraints
  for i = 1:length(dI)
    CIi = CI(i,:); dIi = dI(i);
    c = CIi - (CIi*v)*CE/(CE*v);
    d = dIi - (CIi*v)*dE/(CE*v);
    % normalize the normal vector
    magnitude = sqrt(c*c');
    if (magnitude == 0) 
      if (d < 0)
        error('Infeasible constraint found! Something is wrong!')
      end
    else
      c = c/magnitude;
      d = d/magnitude;
      CIreach = [CIreach; c];
      dIreach = [dIreach; d];
    end
  end
end


% create a linearcon object and return the reachable set
R = linearcon(CEreach,dEreach,CIreach,dIreach);
return