compute_mapping_dha_vr.m

一个matlab的将军模型

function [destination,null_event,time_limit,out_of_bound,terminal] = ...
    compute_mapping_DHA_VR(X0,sttype,Tstamp,srcloc,srccell)

% Compute the `mapping set` from the given polytope in the specified location.
%
% Syntax:
%   "[destination,null_event,time_limit,out_of_bound,terminal] = ...
%			compute_mapping_DHA(X0,sttype,Tstamp,srcloc,srccell,param_file)"
%
% Description:
%
%   The routine compute_mapping_DHA_VR performs the reachability from a initial set X0, finding  the 
%   set of states that can be reached, using variable-rate sampled-data analysis (i.e. discrete-time, difference
%   equations with sample periods that are a function of the continuous-valued state variables).
%   It also finds a conservative approximation of the elapsed time to reach them. 
%
%	 After this, it applies the reset function defined for this transition. As a last
%   operation, it wraps up the set of states with a conservative bounding box, avoiding explosion 
%   in the number of vertices and faces. 
%
%   The inputs to this function are: 
%
%   * "X0": a "linearcon" object representing the initial continuous set
%     for which the mapping is to be computed.
%  
%   * "sttype": type of state for "X0" which is either an `initial` or
%     `face` state (apply the initial displacement or not).
%
%	 * "Tstamp": vector with the time uncertainty to reach the region
% 
%   * "srcloc": source location.
% 
%	 * "srccell": cell region of source
%
%   The outputs of this function are
%
%   % `destination is a cell array with the following fields:
%
%   * dstloc 	: location whose guard was violated during the reachability process
%	 * dstcell 	: Interior region of the dstloc
%   * region 	: linearcon object representing the set of states that violates some guard for a given
%     location.
%   * Tstamp 	: time elapsed interval to reach the region starting at the set X0.
%
%   * "null_event": a boolean flag indicating that the flow pipe computation
%     was terminated because it can be concluded that the subsequent flow
%     pipe segments will remain inside "INV" forever.
%
%   * "time_limit": a boolean flag indicating that the flow pipe computation
%     was terminated because the time limit "max_time" was exceeded.
%
% Implementation:
%   Call the mapping computation routine corresponding to the type of the
%   composite continuous dynamics of all SCSBs for the given location
%   ("fs_lin_map" for `linear` (affine) dynamics, "fs_nonlin_map" for
%   `nonlinear` dynamics, and "clk_map" for `clock` dynamics). %   
%
% See Also:
%   fs_nonlin_map,fs_lin_map,clk_map

global GLOBAL_PIHA
global GLOBAL_APPROX_PARAM
PLOT=1;
% Get the invariant polytope for the given cell location.
INV = return_invariant(srccell);

% Compute the type of the composite dynamics of all SCSBs for the FSM
% state in the source location.

overall_dynamics = check_overall_dynamics(GLOBAL_PIHA.SCSBlocks, ...
    GLOBAL_PIHA.Locations{srcloc}.q);
 
max_time = GLOBAL_APPROX_PARAM.max_time;

null_event = 0;
time_limit = 0;
terminal = {};
out_of_bound = 0

equilibrium=0;
check_equilibrium=0;
check_time_limit=0;

% Get information about analysis region
for k=1:GLOBAL_PIHA.NAR
   c_AR(k,:) = GLOBAL_PIHA.Hyperplanes{k}.c;
   d_AR(k,:) = GLOBAL_PIHA.Hyperplanes{k}.d;
end
AR = linearcon([],[],c_AR,d_AR);


%For now, we treat all systems as nonlinear.
overall_dynamics = 'nonlinear'

switch overall_dynamics,
   
case 'linear',
   
        
case 'clock',
    v = overall_system_clock(GLOBAL_PIHA.SCSBlocks, ...
        GLOBAL_PIHA.Locations{srcloc}.q);
    
    
 case 'nonlinear',
    ode_param = GLOBAL_PIHA.Locations{srcloc}.q;
    check_time_limit = 1;
    
    
 otherwise,
    error(['Unknown dynamics type ''' overall_dynamics '''.'])
    
 end



fprintf(1,'Computing flow pipe segments:')
counter = 0;

global_var;
%PLOT = 0;

switch sttype
   
	case 'init', % initial state
      initial = 1;
      
	case 'face', % face state (a patch)
      initial=0;
          
   otherwise,
      error(['Unknown state type ''' sttype '''.'])
      return
end



%equilibrium = 0; time_limit = 0;Xsim = vertices;
%first = 1; t_total = 0; stop = 0;init_T=T;

% Compute the transformation of the initial segment

time_pointer = 0;
mapping_counter =1;
dst_struc = {};

% Currently, CheckMate only supports one variable-rate clock structure.
for i=1:1
   Pk_reach = X0;
   END_SEARCH = 0;       
	while ~END_SEARCH
      transition_found=0;
      
      % ----- Advance one step of clock  ---------% 
      
      
      if ~isempty(Pk_reach)
         switch overall_dynamics
         
	      	case 'linear'
                
                
            
		    case 'clock'
	   	      
	      
            case 'nonlinear'
               
   	         sample_points_X0 = vertices(Pk_reach);
      	      Pk = linearcon(Pk_reach);          
         	   [Pk,sample_points] = seg_approx_VRClock('overall_system_ode',ode_param,...
            	      Pk,sample_points_X0);
					Pk_reach = Pk;                 
               
          end 	%switch%
          if PLOT        
            drawnow;plot(polyhedron(Pk_reach));  
          end
       
      end	%if%
      
      if ~isempty(Pk_reach)
         
         % -- test if Pk_reach left cell -- %
      	if ~isempty(minus(Pk_reach,return_invariant(srccell)))
               
         	% -- Cycle through possible transitions -- %             
            for m=1:length(GLOBAL_PIHA.Locations{srcloc}.transitions)
                
          		% -- Cycle through each guard cell for each possible transition -- %             
                for k=1:length(GLOBAL_PIHA.Locations{srcloc}.transitions{m}.guard)
                   
                  % -- Test if guard cell was reached by Pk_reach -- %             
						if ~isempty(and(Pk_reach,return_invariant(GLOBAL_PIHA.Locations{srcloc}.transitions{m}.guard(k))))
                     
                     % -- compute intersection between Pk_reach and the guard cell -- %             
                     mapping  = and(Pk_reach,return_invariant(GLOBAL_PIHA.Locations...
                        {srcloc}.transitions{m}.guard(k)));
                     
                     % -- Test if transitions refers to clock used in this loop -- %             
                     if (GLOBAL_PIHA.Locations{srcloc}.transitions{m}.clock == i)
                        
                        % -- Compute location number of the destination (empty if terminal) -- % 
                            destination_sfm				= GLOBAL_PIHA.Locations{srcloc}.transitions{m}.destination;
						    [dstloc,destin_location]	= find_location(srcloc,destination_sfm,m);
                        
                        if isempty(dstloc)
                           a=1;
                        end
                      
                        % -- Test if location is terminal -- %             
                        if ~isempty(dstloc)
                           
	                     % -- Apply reset to the reachable set-- %             
                           mapping_pos_reset 	= apply_reset(mapping,GLOBAL_PIHA.Locations{srcloc}.transitions{m}.destination...
                              ,GLOBAL_PIHA.Locations{srcloc}.transitions{m});
                           if PLOT
                               plot(mapping_pos_reset,[1 0 0]);
                           end
                        
		                     % -- cycle through all the interior cells for the location of the destination -- %             
                        	for d=1:length(GLOBAL_PIHA.Locations{dstloc}.interior_cells)
                              
                              cell_number = GLOBAL_PIHA.Locations{dstloc}.interior_cells(d);
                              
				                  % -- Test if mapping intersect given cell of the destination -- %             
                              if ~isempty(and(mapping_pos_reset,return_invariant(cell_number)))
                                 
                               	dst_struc{mapping_counter}.dstloc  	= dstloc;
	                           	dst_struc{mapping_counter}.mapping	= and(mapping_pos_reset,return_invariant(cell_number));
   	                        	dst_struc{mapping_counter}.dstcell 	= cell_number;
      	                  		dst_struc{mapping_counter}.Tstamp   = [inf_pointer(i) sup_pointer(i)];
                                dst_struc{mapping_counter}.transition_theta=0;
                                 if ~strcmp(overall_dynamics,'clock')
                                     %Only perform bounding box if there is more than one mapping
                                     %('if then' added by JPK 1/30/02)
                                     if length(dst_struc{mapping_counter}.mapping)>1
                                        mapping_pos_bb	= bounding_SD_box(dst_struc{mapping_counter},param_file);   
                                    else
                                        mapping_pos_bb	=   dst_struc{mapping_counter}.mapping;
                                    end    
                                else
                                    mapping_pos_bb	=   dst_struc{mapping_counter}.mapping;   
                                 end
            	               	dst_struc{mapping_counter}.mapping	= mapping_pos_bb;
               	            	mapping_counter = mapping_counter + 1;
                  	         end %if%
                           end %for%
                        else
                           %Else add to terminal list if the terminal state has not already been added.
                           found_terminal=0;
                           for q=1:length(terminal)
                               if terminal{q}==destination_sfm
                                   found_terminal=1;
                               end
                           end   
                           if ~found_terminal
                                terminal{length(terminal)+1} = destination_sfm;
                           end
                           
                        end %if%
                  	end%if%
                  	transition_found =1;
                  	temp = minus(Pk_reach,mapping);
          				if ~isempty(temp)
            				Pk_reach  = temp{1};
         				else
            				Pk_reach = temp;
         				end
                       
						end	%if%
		     		end	%for%        
				end	%for%
			end	%if%
      end	%if%
      
      % If guard violated, collect information about the cells in the source location
      if transition_found |(initial & ~isempty(Pk_reach))
         
         % If reachable set is not empty