hyperrectangle.m

CheckMate is a MATLAB-based tool for modeling, simulating and investigating properties of hybrid dyn

function SEG=hyperrectangle(sys_eq,ode_param,X0,t0,tf,Pcon,Siminf,SP0,SPf)


% Compute convex hull from these extreme points
% CH = polyhedron(SP0 | SPf);
[CE,dE,C,d] = linearcon_data(linearcon(polyhedron(SP0 | SPf)));

global GLOBAL_APPROX_PARAM

[CE,dE,CI,dI] = linearcon_data(X0);
[CPE,dPE,CPI,dPI] = linearcon_data(Pcon);

%This section has been altered (3/2002) to accomodate lower dimensional linearcon
%objects.  This was introduced specifically for parametric analysis.
% (1) inequality constraints

    [CIn,CIm]=size(CI);
    [CPIn,CPIm]=size(CPI);
    xdim=dim(X0);
    pdim=dim(Pcon);

    % Remember the dimensions of the state space.
    dimension=xdim;
    fminCI=[zeros(1,xdim) 1 zeros(1,pdim); zeros(1,xdim) -1 zeros(1,pdim);];%CIm+1
    fmindI=[tf;-t0];
    
    if ~isempty(CI)
        fminCI=[fminCI;CI,zeros(CIn,pdim+1)];%CPIm+1
        fmindI=[fmindI;dI];
    end;
    
    if ~isempty(CPI)
        fminCI=[fminCI;zeros(CPIn,xdim+1),CPI];%CIm+1
        fmindI=[fmindI;dPI];
    end;

    % (2) equality constraints
    [CEn,CEm]=size(CE);
    [CPEn,CPEm]=size(CPE);

    fminCE=[];
    fmindE=[];
    if ~isempty(CE)
        fminCE=[CE,zeros(CEn,pdim+1)];%CPIm+1
        fmindE=dE;
    end;
    if ~isempty(CPE) 
        if ~isempty(fminCE)
            fminCE=[fminCE;zeros(CPEn,xdim+1) CPE];%CIm+1
            fmindE=[fmindE;dPE];
        else
            fminCE=[zeros(CPEn,xdim+1) CPE];%CIm+1
            fmindE=[dPE];
        end;
    end;

% (3) Set the optimization parameters

tolerance = GLOBAL_APPROX_PARAM.func_tol;
maxiter = GLOBAL_APPROX_PARAM.max_iter;
if (tf-t0)>0
    actual_tolerance = (tf-t0)*tolerance;
else
    actual_tolerance = tolerance;
end


options = optimset;
options.TolFun =actual_tolerance;%f(3) of foptions
options.MaxIter = maxiter*size(C,2);%f(14) of foptions
options.Display=['off'];
options.Diagnostics=['off'];
options.LargeScale=['off'];


if GLOBAL_APPROX_PARAM.verbosity>=1
   fprintf('\n\nCompute flowpipe segment with %i faces:\n',size(C,1))
end;

% (4) Optimize for each direction

for l = 1:size(C,1)
  n_vector = C(l,:)';
  dlmax = Inf;
  
    % compute the best initial state 
    % Store in " tindex"  the index of the time elements of  the vector Siminf(1).T that are greater than t0
    %
    [ts,tindex]=setdiff(Siminf(1).T.*(Siminf(1).T>=t0),0);
    % Computer the max arg i C'x(i)
    [minnx,nx_index]=min(-n_vector'*Siminf(1).X(tindex,:)');
    Pinit=Siminf(1).P0;
    Xinit=Siminf(1).X0;
    Tinit=Siminf(1).T(tindex(nx_index));
    for k0=2:length(Siminf)
        [ts,tindex]=setdiff(Siminf(k0).T.*(Siminf(k0).T>=t0),0);
        [nx,nx_index]=min(-n_vector'*Siminf(k0).X(tindex,:)');
        if nx<minnx
            minnx=nx;
            Pinit=Siminf(k0).P0;
            Xinit=Siminf(k0).X0;
            Tinit=Siminf(k0).T(tindex(nx_index));
        end;
    end;
    if GLOBAL_APPROX_PARAM.verbosity==1
        fprintf('%3i ',l)
        if mod(l,20)==0
            fprintf('\n')
        end
    end;
    if GLOBAL_APPROX_PARAM.verbosity>=2
        fprintf('\n%4i:',l);
        fprintf('x0=[');
        fprintf(' %4.2f',Xinit);
        if isempty(Pinit)
            fprintf(']')
        else
            fprintf(']\tp=[');
            fprintf(' %4.2f',Pinit);
            fprintf(']');
        end;
    end;
 
    
    [Xopt,dl] = fmincon('stretch_func_ode',[Xinit;Tinit;Pinit], fminCI,fmindI, fminCE,fmindE,[],[],[],options,...
             sys_eq,ode_param,n_vector,t0,tf,dimension);
    

    d(l,1)=-dl;
end

% >>>>>>>>>>>> Create linearcon Object -- DJ -- 06/30/03 <<<<<<<<<<<<
% Added by Dong Jia 
SEG=linearcon([],[],C,d);
% >>>>>>>>>>>> -------------- end (Create linearcon Object) --------------- <<<<<<<<<<<<