iauto_part.m

一个matlab的将军模型

            entry_region{m}.state{l}.polytope = partition{l};
            entry_region{m}.state{l}.mapping = {};
            entry_region{m}.state{l}.children = [];
            entry_region{m}.state{l}.null_event = -1;
            entry_region{m}.state{l}.time_limit = -1;
            entry_region{m}.state{l}.out_of_bound = -1;
            entry_region{m}.state{l}.indeterminate = 0;
            entry_region{m}.state{l}.non_reachable = 1;
            entry_region{m}.state{l}.terminal = {};
            entry_region{m}.state{l}.visited  = 0;
            entry_region{m}.state{l}.split = 0;
            % >>>>>>>>>>>> Adding Field destination (states) -- DJ -- 06/30/03 <<<<<<<<<<<<
            % The field destination is add to record the reachability results.
            % Added by Dong Jia
            entry_region{m}.state{l}.destination = {};
            % >>>>>>>>>>>> -------------- end (Adding Field destination (states)) --------------- <<<<<<<<<<<<
        end 
        % >>>>>>>>>>>> -------------- end (Initial Partition) --------------- <<<<<<<<<<<<
    end 
    if isempty(entry_region)
        GLOBAL_AUTOMATON{k} = [];
    else
        GLOBAL_AUTOMATON{k}.interior_region = entry_region;
    end 
end 
fprintf(1,'\n')
if init_reach
    %% Refine partition for invariant faces in each location further using
    %% reachability analysis.
    for m = 01:length(GLOBAL_PIHA.InitialConditions)
        % Compute sequence of reachable polytopes for each initial state
        init_loc = GLOBAL_PIHA.InitialConditions{m}.initialLocation;
        init_cell=GLOBAL_PIHA.InitialConditions{m}.initialCells;
        for n=1:length(GLOBAL_AUTOMATON{init_loc}.initstate)
            fprintf(1,'Reachability analysis for loc %d , Cell region %d initial state %d:\n',init_loc,init_cell,n)
            [reachable,time_limit] = ...
                flow_reach(GLOBAL_AUTOMATON{init_loc}.initstate{n}.polytope,init_loc,init_cell);
            % Intersect (subtract) each reachable polytope with (from) all states
            % on the same face.
            % Skip the first state in the reachable set which is always an
            % initial state.
            GLOBAL_AUTOMATON{init_loc}.initstate{n}.time_limit = time_limit;
            GLOBAL_AUTOMATON{init_loc}.initstate{n}.mapping = reachable(1).mapping;
            % >>>>>>>>>>>> Recording Reachability Results (initial state) -- DJ -- 06/30/03 <<<<<<<<<<<<
            % Recording the reachability results for the initial state.
            % Added by Dong Jia to record complet information from
            % reachability analysis
            GLOBAL_AUTOMATON{init_loc}.initstate{n}.destination=reachable(1).destination;
            GLOBAL_AUTOMATON{init_loc}.initstate{n}.null_event=reachable(1).null_event;
            GLOBAL_AUTOMATON{init_loc}.initstate{n}.time_limit=reachable(1).time_limit;
            GLOBAL_AUTOMATON{init_loc}.initstate{n}.out_of_bound=reachable(1).out_of_bound;
            GLOBAL_AUTOMATON{init_loc}.initstate{n}.terminal=reachable(1).terminal;
            % >>>>>>>>>>>> -------------- end (Recording Reachability Results (initial state)) --------------- <<<<<<<<<<<<
            
            fprintf(1,'\n Finished preliminary reachability analysis. \n');
            % save reachable_data reachable;
            prel_refine(reachable);
        end %for n
    end % for m
end %if init_reach

%Section ends here ---------------------------------------------------------------------

return


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

function prel_refine(reachable)


global GLOBAL_AUTOMATON;
global GLOBAL_PIHA;

fprintf(1,'\n Starting preliminary refinement.\n');
len_reachable = length(reachable);

%Loop through the reachable set
for m = 2:length(reachable)
    rlocation = reachable(m).location;
    [dum1,dum2,int_reg] = intersect(reachable(m).cell,GLOBAL_PIHA.Locations{rlocation}.interior_cells);
    polytope = reachable(m).polytope;
    
    %Get current states for this face
    state = GLOBAL_AUTOMATON{rlocation}.interior_region{int_reg}.state;
    refinement = {};
    
    %Loop through current states and break apart based on reachable
    %set mapping.
    for n = 1:length(state)
        temp = state{n};
        overlap = polytope & state{n}.polytope;
        if ~isempty(overlap)
            [CE,dE,CI,dI]=linearcon_data(overlap);
            box = grow_polytope_for_iautopart(CE,dE,CI,dI,rlocation,state{n}.polytope);
            temp.polytope = (box&return_invariant(reachable(m).cell));
            temp.mapping  = reachable(m).mapping;
            temp.non_reachable = 0;
            % >>>>>>>>>>>> Recording Reachability Results (Intermedia States) -- DJ -- 06/30/03 <<<<<<<<<<<<
            % Recording reachability results for intermedia results.
            % Added by Dong Jia to record all information from reachability
            % analysis.
            temp.destination=reachable(m).destination;
            temp.null_event=reachable(m).null_event;
            temp.time_limit=reachable(m).time_limit;
            temp.out_of_bound=reachable(m).out_of_bound;
            temp.terminal=reachable(m).terminal;
            % >>>>>>>>>>>> -------------- end (Recording Reachability Results (Intermedia States)) --------------- <<<<<<<<<<<<
            % >>>>>>>>>>>> Split Flag -- DJ -- 06/30/03 <<<<<<<<<<<<
            % Add a flag split to indicate whether a state is created by splitting a
            % previous state. Added by Dong Jia
            temp.split = 0;
            % >>>>>>>>>>>> -------------- end (Split Flag) --------------- <<<<<<<<<<<<
            refinement{length(refinement)+1} = temp;
            fprintf(1,'\n  state %d/%d created.\n',m,len_reachable);
        end
    end
    
    if length(refinement)>1
        for n=1:length(refinement)
            refinement{n}.split=1;
        end
    end
    
    % This part find new reachable regions and add them to the structure called 'refinement'. 
    % At the end, the structure 'refinement' is copied to 
    % GLOBAL_AUTOMATON{LOCATION}.interior_Region{REGION}.state
    % The cycle then repeats itself for length of structure 'reachable'.
    % >>>>>>>>>>>> Getting States for Current Cell -- DJ -- 06/30/03 <<<<<<<<<<<<
    % Get the state list for current cell to refine.
    % Added by Dong Jia
    curr_states = GLOBAL_AUTOMATON{rlocation}.interior_region{int_reg}.state;
    % >>>>>>>>>>>> -------------- end (Getting States for Current Cell) --------------- <<<<<<<<<<<<
    % >>>>>>>>>>>> State Refinement (while part)-- DJ -- 06/30/03 <<<<<<<<<<<<
    % In the while loop, all states will be checked to see if they intersect
    % with each other. If so, the two intersected states will be partitioned.
    % The following while-end loop is added by Dong Jia
    k=0;
    while 1
        k=k+1;
        if k>length(refinement)
            break;
        end
        % >>>>>>>>>>>> -------------- end (State Refinement (while part)) --------------- <<<<<<<<<<<<
        % >>>>>>>>>>>> Eliminating non_reachable-- DJ -- 06/30/03 <<<<<<<<<<<<
        % The field non_reachable is not used to judge whether one state might need to be
        % refined. The if-end is commented by Dong Jia
        %       if ~(GLOBAL_AUTOMATON{rlocation}.interior_region{int_reg}.state{j}.non_reachable)
        % >>>>>>>>>>>> -------------- end (Eliminating non_reachable) --------------- <<<<<<<<<<<<
        intersect_flag=1;
        for j=1:length(curr_states)
            poly=curr_states{j}.polytope;
            if isfeasible(poly,refinement{k}.polytope)
                inter=refinement{k}.polytope&poly;
                intersect_flag=0;
                % The if-elseif-else-end block 
                % is to consider special cases when one set is a subset of the other one.
                if issubset(refinement{k}.polytope,poly)
                    %   poly contain a state in the refinement list. in
                    %   this case poly will be partitioned. All states
                    %   generated from poly-inter are appended to the
                    %   curr_states list and curr_state{j} is updated by
                    %   inter.
                    GA_not=poly-inter;
                    for i=1:length(GA_not)
                        curr_states{length(curr_states)+1}=curr_states{j};
                        curr_states{length(curr_states)}.polytope=GA_not{i};
                        % >>>>>>>>>>>> Updating split (Case 1) -- DJ -- 06/30/03 <<<<<<<<<<<<
                        % Updating the field split for newly created states
                        % Added by Dong Jia
                        if ~curr_states{j}.non_reachable
                            curr_states{length(curr_states)}.split=1;
                        end
                        % >>>>>>>>>>>> -------------- end (Updating split (Case 1)) --------------- <<<<<<<<<<<<
                    end
                    % >>>>>>>>>>>> Update Information for inter (case 1)-- DJ -- 06/30/03 <<<<<<<<<<<<
                    % Updating reachability information for inter. If curr_state{j} is not a
                    % reachable set computed from  a initial set, update the reachability
                    % infor from refinement{k}. Otherwise, update the reachability infor by the
                    % intersection of the corresponding fields in curr_state{j} and
                    % refinement{k}.
                    % Added by Dong Jia
                    if curr_states{j}.non_reachable
                        curr_states{j}=refinement{k};
                        curr_states{j}.polytope=inter;
                    else
                        curr_states{j}.polytope=inter;
                        if ~isempty(refinement{k}.mapping)
                            if ~isempty(curr_states{j}.mapping)
                                tmp_mapping={};
                                u=1;
                                for s=1:length(curr_states{j}.mapping)
                                    for t=1:length(refinement{k}.mapping)
                                        if isfeasible(refinement{k}.mapping{t},curr_states{j}.mapping{s});
                                            u=u+1;
                                            tmp_mapping{u}=refinement{k}.mapping{t} & curr_states{j}.mapping{s};
                                        end
                                    end
                                end
                                curr_states{j}.mapping = tmp_mapping;
                                tmp_dest={};
                                u=0;
                                for s=1:length(curr_states{j}.destination)
                                    for t=1:length(refinement{k}.destination)
                                        if (refinement{k}.destination{t}.location==curr_states{j}.destination{s}.location) && ...
                                                (refinement{k}.destination{t}.cell==curr_states{j}.destination{s}.cell) && ...
                                                isfeasible(refinement{k}.destination{t}.mapping{1},curr_states{j}.destination{s}.mapping{1})
                                            u=u+1;
                                            tmp_dest{u} = curr_states{j}.destination{s};
                                            tmp_dest{u}.mapping{1}= refinement{k}.destination{t}.mapping{1}&curr_states{j}.destination{s}.mapping{1};
                                            break;
                                        end
                                    end
                                end
                                curr_states{j}.destination=tmp_dest;
                            else
                                curr_states{j}.mapping = refinement{k}.mapping;
                                curr_states{j}.destination=refinement{k}.destination;
                            end
                        end
                        curr_states{j}.null_event=curr_states{j}.null_event & refinement{k}.null_event;
                        curr_states{j}.time_limit=curr_states{j}.time_limit & refinement{k}.time_limit;
                        curr_states{j}.out_of_bound=curr_states{j}.out_of_bound & refinement{k}.out_of_bound;
                        tmp_terminal = {};
                        u=0;
                        for s=1:length(curr_states{j}.terminal)
                            for t=1:length(refinement{k}.terminal)
                                if curr_states{j}.terminal{s}==refinement{k}.terminal{t}
                                    u=u+1;
                                    tmp_terminal{u}=curr_states{j}.terminal{s};
                                    break;
                                end
                            end
                        end
                        curr_states{j}.terminal=tmp_terminal;
                        curr_states{j}.split=refinement{k}.split;
                    end
                    % >>>>>>>>>>>> -------------- end (Update Information for inter (case 1)) --------------- <<<<<<<<<<<<
                    break;
                end
            end %if isfeasible...
            
        end %for j=1....