image thumbnail

pmap - Parameter Space Stability Mapping Suite

by

 

30 Mar 2011 (Updated )

Finds the Hurwitz stability region in a parameter space of a continuous time system.

pmap_generate(spoly, varargin)
function Test = pmap_generate(spoly, varargin)
% Test = pmap_generate(spoly,varargin)
% 
% Version 1.0
% 3-18-11
%
% Generates the structure needed to map the parameter space for stability
% of dynamical systems. The mapping alogrithm uses Sign-Definite
% Decomposition and the contstruction of vertex polynomials reminiscent of
% Kharitonov's theorem. For a detailed look at the approach used here, see 
%
% [1] L.H. Keel and S.P. Bhattacharyya, "Fixed order multivariable 
%     controller synthesis: A new algorithm," Proceedings of the 47th IEEE 
%     Conference on Decision and Control, Cancun, Mexico, December, Cancun,
%     Mexico: 2008.
%
% [2] M.J. Knap, L.H. Keel, and S.P. Bhattacharyya, "Robust Stability of
%     Complex Systems with Applications to Performance Attainment
%     Problems," Proceedings of the 2010 American Control Conference,
%     Baltimore, MD, USA: 2010, pp. 3907-3912
%
%
% Copyright (C) 2011 Michael Knap

% FIXME: Throroughly introduce the use of the different parameter value pairs

%% Parse
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
p=inputParser;
p.addRequired('spoly',@validatepoly);  % required

p.addParamValue('title',[],@ischar);
p.addParamValue('threshold',16,@(x) x>4 && mod(x,1)==0 );
p.addParamValue('debug',3,@(d) isscalar(d) && mod(d,1)==0 );
p.addParamValue('phasemargin',[0 0], @(p) all(size(p)==[1,2]) );
p.addParamValue('gainmargin', [1 1], @(k) all(size(k)==[1 2]) );
p.addParamValue('testmethod','vertex_compan',...  % new default does all of them
    @(str) any(strcmpi(str,...
    {'vertex_compan','phase2','discrete','vertex','gain','phase','gain+phase','none'})));
    % The different types of tests.
    % The vertex test is complete, it was the first method
    % the vertex_compan test is about 3 times faster, but there should be
    % some more tests for the condition number of the companion matrices
    % The others are left-overs, and probably will be removed soon
    
p.CaseSensitive=true;
p.KeepUnmatched=true;
p.FunctionName='pmap_generate';
p.parse(spoly,varargin{:})

Test=p.Results;
if Test.debug>3
    disp('-------------------------------------------------------------------')
	disp('Using pmap_generate')
	disp('Parse Results (p.Results) : ')
	disp(p.Results)
	disp('Unmatched Parse Results (p.Unmatched) :')
	disp(p.Unmatched)
	disp('-------------------------------------------------------------------')
end
tgen=tic;
fprintf('Generating the test structure...\n');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Parsed and validated, get our list of parameters from the symbolic variables
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Here, we are removing the specific symbolic variables 's', 'K', and 'm' 
% from the list of symbolic variables. The remaining variables must be our 
% parameters for the given system. 

% This used to be done with a for loop, but I tried to hack it the 
% "matlab vectorized" way. 
svars=arrayfun(@char,symvar(Test.spoly),'UniformOutput',false);
indices=cellfun(@strcmp,{'s','m','K'},repmat({svars},1,3),'UniformOutput',false);
svars(logical(sum(vertcat(indices{:})) ))=[];

% Old loop code for giggles: 
% svars=symvar(spoly);
% sVars=cell(1,length(svars));
% for k=1:length(svars)
%     sVars{k}=char(svars(k));
%     if any(strcmp(char(svars(k)),{'s','K','m'})) % these are not parameters
%         svars(k)=NaN;
%     end
% end
Test.dims=svars;          % cell array of strings
Test.params=sym(svars);   % array of symbolic variables
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Now we need to get and setup our parameter space:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FIXME: I used to request the parameter space it it were not given, now 
% I am assigning a default space. The potential problem with this is the 
% creation of unneeded polynomials from unused orthants. For now, I am 
% going to leave it this way. 
ndims=length(Test.dims);
pts=repmat([-16;16],1,ndims);  % powers or 2 are good defaults for bisection

for k=1:ndims
    if isfield(p.Unmatched,Test.dims{k})
        disp(['Setting ' Test.dims{k} ' range : [' num2str(p.Unmatched.(Test.dims{k})) '].']);
        pts(:,k)=p.Unmatched.(Test.dims{k})';
    end
end
Test.pts=sort(pts);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Now we construct our orthants and assumptions:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Test=GetOrthants(Test);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% For the phase margin intervals...
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% These could possibly be eliminated if they were hard-coded in the
% actual test function. For now I want to leave them here. 
Test.intervals=[ 1 2 3 4];
Test.intervalassums={
        {'Type::Negative' 'Type::Negative'}   % interval I
		{'Type::Negative' 'Type::Positive'}    % interval II
		{'Type::Positive' 'Type::Positive'}    % interval II
		{'Type::Positive' 'Type::Negative'}    % interval III
		};
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Create a function handle that will can evaluate as a numeric polynomial
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This is exteremely useful in testing, debugging, and verifying results.
% With this we can evaluate specific points in the parameter space and
% margin space 
Test.polyfcn=matlabFunction(fliplr(coeffs(Test.spoly,sym('s'))).',...
    'vars', [Test.dims { 'm' 'K' }] ); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Create the vertex polynomials (and companion matrices for new version)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This is the real meat of the algorithm. There are lots of neccessary
% sub-functions called from here. 
Test=CreateVertexPolynomials(Test);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Tidy Up... create a title for this test if not provided
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% We also want to provide a timestamp and version stamp
if Test.debug>4
    ttitle=tic;
    fprintf('     Stamping and setting title...');
end
Test.version.generate=1.0;
Test.date.generate=datestr(now);

if exist('linewrap.m','file')==2
    printpoly=char(linewrap(char(collect(Test.spoly,'s')),80));
else
    if Test.debug > 3
        warning('pmap_generate:title','   Cannot do linewrap on polynomial')
    end
    printpoly=char(collect(Test.spoly,'s'));
end
if isempty(Test.title)
    Test.title={ 'Hurwitz Stability Set of ';
		printpoly };
        % The title is better with just the polynomial. 
        % The phase and gain margins are printed with legend entries now.
        
		%['Gain margin  : [' num2str(Test.gainmargin) ']']
		%['Phase margin : [' num2str(Test.phasemargin) ']']
		%}; 
end
if Test.debug>4
    fprintf('(%4.2f s)\n',toc(ttitle));
end
Test.t.generate=toc(tgen);
fprintf('Done generating test structure (%4.2f s)\n',Test.t.generate);
if Test.debug>4
    disp(Test);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Supporting Functions %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% bool = validatepoly(spoly)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% this is called from the input parser to validate the input polynomial
function bool=validatepoly(spoly)
assert(isa(spoly,'sym'),'bkkGenerate2:error','The polynomial must be symbolic.');

svars=arrayfun(@char,symvar(spoly),'UniformOutput',false);
% vars=symvar(spoly);
% for k=1:length(vars)
%     Vars{k}=char(vars(k));
% end

assert(any(strcmp('s',svars)), 'bkkGenerate2:error','Polynomial must be a polynomial in ''s''.');

bool=true;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Test = GetOrthants(Test)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This has been rewritten with parallel code
% Rewrite 3-3-11 - passing and returning whole test structure. 
% Rewrite 3-18-11 ... adding support for phase intervals in the dimension
% of the orthants

% FIXME : This can be made MUCH cleaner

function Test = GetOrthants(Test)
torth=tic;
if Test.debug>=1
    fprintf('     Getting orthants and assumptions... ');
end
Points=Test.pts;
dimsToDivide= 0 < Points(2,:) & 0 > Points(1,:);
numOrths=2^sum(dimsToDivide);
if numOrths>1
    orthants=cell(1,numOrths);   % the row dimentions of orthants will be the phase intervals
    
    % Constructs the upper and lower parts of each dim to divide
    lower=cell(1,length(dimsToDivide));
    upper=cell(1,length(dimsToDivide));
    
    for k=1:length(dimsToDivide)
        if dimsToDivide(k)
            lower{k}=[Points(1,k) 0]';
            upper{k}=[0 Points(2,k)]'; 
        end
    end
    
    % uses binary to construct orthants from upper and lower parts
    for k=1:numOrths
        b=bitget(k-1,length(dimsToDivide):-1:1);
        for jj=1:length(dimsToDivide)
            if dimsToDivide(jj)==1
                if b(jj)==0
                    orthants{1,k}(:,jj)=lower{jj};
                else
                    orthants{1,k}(:,jj)=upper{jj};
                end
            else
                orthants{1,k}(:,jj)=Points(:,jj);
            end
        end
    end
    
    % Now we define the assumptions for each orthant
    assums=cell(1,numOrths);
    ass=cell(1,length(dimsToDivide));
    for k=1:numOrths
        for jj=1:length(dimsToDivide)
            tmp=min(orthants{:,k}(:,jj));
            if tmp<0
                ass{jj}='Type::Negative';
            else
                ass{jj}='Type::Positive';
            end
        end
        assums{k}=ass;
    end
else
    orthants{:,1}=Points;
    ass=cell(1,length(dimsToDivide));
    for jj=1:length(Points)
        tmp=min(Points(:,jj));
        if tmp<0
            ass{jj}='Type::Negative';
        else
            ass{jj}='Type::Positive';
        end
    end
    assums{1}=ass;
end
Test.orthants=orthants;
Test.assums=assums;
Test.t.orth=toc(torth);
if Test.debug >=1
    fprintf(' %d orthants found (%4.2f s).\n',length(orthants),Test.t.orth)
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% FIXME from here down
%% Test = CreateVertexPolynomials(Test)  (Pcompan version)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This is the real meat of the algorithm with my addition of using the 
% orthant-based approach. The use of orthant-based assumptions has proven
% to be quite a lot of work, but the resulting sets are much sharper.

% One can see a difference between the orthant based approach and the
% translation approach by testing the same polynomial, one with translation
% and one without. The translation makes the polynomials much more complex,
% hence the roots are more susceptible to the distance of the vertex
% polynmomails and the parameter space require extended bisection. 

% I have also constructed the orthant-based approach for the phase-margin. 

% I have also moved to the use of evaluatable companion matrices. This has
% sped up the calculations by about a factor of 3 for the systems I have tested. 

function Test = CreateVertexPolynomials(Test)
tvert=tic;
if Test.debug>0
    fprintf('     Creating per-orthant vertex polynomials...');
end
parfor ii=1:4   % ii runs over the different phase intervals
    tinterval=tic;
    if Test.debug>1
        fprintf('\n          Phase interval %d...',ii);
    end
    Pcompan{ii}=PerInterval(Test,ii);
    %Test.P{ii}=P;
    if Test.debug>1
        fprintf('created %d polynomials (%4.2f s).',...
            8*length(Test.orthants),toc(tinterval));
    end
end
Test.Pcompan=vertcat(Pcompan{:}); % the new syntax should be Pcompan{interval, orthant};
Test.t.vert=toc(tvert);
if Test.debug>0
    fprintf('(%4.2f s)',Test.t.vert);
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Pcompan=PerInterval(Test)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Trying to paralellize this damn code...
function Pcompan=PerInterval(Test,ii)
Pcompan=cell(1,length(Test.orthants));

for k=1:length(Test.orthants)
    % FIXME: this is a carry-over from a previous version.
    % Can I write this differently?
    
    params=Test.params;
    assum=Test.assums{k};
    intassum=Test.intervalassums{ii};
    %% First, we get our symbolics with the proper assumptions
    for j=1:length(params)
        feval(symengine,'assume',params(j),assum{j});
    end
    syms s K m real; % changed Theta to m 3-3-11
    %% Real and Imaginary Decomposition
    a=feval(symengine,'Re',Test.spoly);
    b=feval(symengine,'Im',Test.spoly);
    %% Even and Odd Decomposition
    [aeven aodd]=EvenOddTerms(a);
    [beven bodd]=EvenOddTerms(b);
    %% Positive and Negative Descomposition
    %keyboard
    [aevenpos aevenneg]=GetPosNeg(aeven,params,assum,intassum);
    [aoddpos aoddneg]=GetPosNeg(aodd, params,assum,intassum);
    [bevenpos bevenneg]=GetPosNeg(beven, params,assum,intassum);
    [boddpos boddneg]=GetPosNeg(bodd, params, assum,intassum);
    %% Construct the vertex polynomial pieces
    %keyboard
    [aevenplus aevenminus]=DoEvals(aevenpos,aevenneg,params,assum,Test.intervals(ii));
    [aoddplus aoddminus]=DoEvals(aoddpos, aoddneg,params,assum,Test.intervals(ii));
    [bevenplus bevenminus]=DoEvals(bevenpos, bevenneg, params, assum,Test.intervals(ii));
    [boddplus boddminus]=DoEvals(boddpos, boddneg,params,assum,Test.intervals(ii));
    %% Put the pieces together
    p1=aevenminus+1i*boddplus+1i*(bevenminus-1i*aoddminus);
    p2=aevenminus+1i*boddplus+1i*(bevenplus-1i*aoddplus);
    p3=aevenplus+1i*boddminus+1i*(bevenplus-1i*aoddplus);
    p4=aevenplus+1i*boddminus+1i*(bevenminus-1i*aoddminus);
    p5=aevenminus+1i*boddminus+1i*(bevenminus-1i*aoddplus);
    p6=aevenminus+1i*boddminus+1i*(bevenplus-1i*aoddminus);
    p7=aevenplus+1i*boddplus+1i*(bevenplus-1i*aoddminus);
    p8=aevenplus+1i*boddplus+1i*(bevenminus-1i*aoddplus);
    
    %% Create Functions:
    % Create our extreme parameters and gain
    minparams=cell(1,length(params));
    maxparams=cell(1,length(params));
    for i=1:length(params)
        minparams{i}=['min' char(params(i))];
        maxparams{i}=['max' char(params(i))];
    end
    
    %% Return our set:
    switch lower(Test.testmethod)
        
        case 'test1'
            % was doing some testing with this switch
        
        case 'vertex'
            % this was the original way, but the calc code had to loop
            % through each polynomial to find roots
            Pcompan{k}={...
                matlabFunction(fliplr(coeffs(p1,'s')).',...
                'vars',[minparams maxparams {'minm' 'maxm' 'minK' 'maxK'}]),...
                matlabFunction(fliplr(coeffs(p2,'s')).',...
                'vars',[minparams maxparams {'minm' 'maxm' 'minK' 'maxK'}]),...
                matlabFunction(fliplr(coeffs(p3,'s')).',...
                'vars',[minparams maxparams {'minm' 'maxm' 'minK' 'maxK'}]),...
                matlabFunction(fliplr(coeffs(p4,'s')).',...
                'vars',[minparams maxparams {'minm' 'maxm' 'minK' 'maxK'}]),...
                matlabFunction(fliplr(coeffs(p5,'s')).',...
                'vars',[minparams maxparams {'minm' 'maxm' 'minK' 'maxK'}]),...
                matlabFunction(fliplr(coeffs(p6,'s')).',...
                'vars',[minparams maxparams {'minm' 'maxm' 'minK' 'maxK'}]),...
                matlabFunction(fliplr(coeffs(p7,'s')).',...
                'vars',[minparams maxparams {'minm' 'maxm' 'minK' 'maxK'}]),...
                matlabFunction(fliplr(coeffs(p8,'s')).',...
                'vars',[minparams maxparams {'minm' 'maxm' 'minK' 'maxK'}])};
            
            
		otherwise
			% now, the calc code can calculate all the roots of the 8
			% polynomials simultaneously with this block diagonal companion
			% matrix
			Pcompan{k}=matlabFunction(blkdiag(...
                feval(symengine,'linalg::companion',p1,'s'),...
                feval(symengine,'linalg::companion',p2,'s'),...
                feval(symengine,'linalg::companion',p3,'s'),...
                feval(symengine,'linalg::companion',p4,'s'),...
                feval(symengine,'linalg::companion',p5,'s'),...
                feval(symengine,'linalg::companion',p6,'s'),...
                feval(symengine,'linalg::companion',p7,'s'),...
                feval(symengine,'linalg::companion',p8,'s')),...
                'vars',[minparams maxparams {'minm' 'maxm' 'minK' 'maxK'}]);
            %         P{k}=vpoly;
    end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Supporting Function EvenOddTerms
function [ev od]=EvenOddTerms(p)
syms s ;
ev=sym(0);
od=sym(0);
degree=feval(symengine,'degree',p,s);
for n=degree:-1:0
	this=feval(symengine,'coeff',p,s,n);
	if mod(n,2)==1  % odd
		od=od + (this)*s^n;
	else
		ev=ev + (this)*s^n;
	end
end

end

%% Supporting Function GetPosNeg
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Adding support to get positivity of cos and sin on certain intervals

function [pos neg]=GetPosNeg(p,params,assums,intervalassums)
		%keyboard;
		%% Make our parametric assumptions :
		syms s m K real;
		syms K positive;
		% Or can we skip the parms and assums ?
		for k=1:length(params)
			feval(symengine,'assume',params(k),assums{k});
		end
		% if we have assigned assumptions for the intervals, 
		%keyboard;
        if length(intervalassums)==2
            feval(symengine,'assume'    ,'sin(m)',intervalassums{1});
            feval(symengine,'assumeAlso','cos(m)',intervalassums{2});
        end
		
		pos=sym(0);
		neg=sym(0);
		%p=simple(p); % added Thu. May 13
		degree=feval(symengine,'degree',p,s);
		% NOTE : this seems to have fixed things, but can we make this
		% NOTE : better, or parallelize it ?
		for n=degree:-1:0
			q=feval(symengine,'coeff',p,s,n);
			
			%% Algorithm with assumptions :
			Nterms=feval(symengine,'nterms',q);
			True=feval(symengine,'TRUE');
			False=feval(symengine,'FALSE');
			%Unknown=feval(symengine,'UNKNOWN');
			
			for i=1:Nterms
				t=feval(symengine,'nthmonomial',q,i);
% 				trigexpr=regexp(char(t),'(sin|cos)(\().*(\))+','match');
% 				if ~isempty(trigexpr)
% 					feval(symengine,'assume',trigexpr,'Type::Positive');
% 				end
				bool=feval(symengine,'is',t,'Type::Positive');
				if bool==True
					pos=pos + t*s^n;
				elseif bool== False
					neg=neg +t*s^n;
                else
                    cprintf('red', '\nCan not get positivity of %s.\nLine 465\n',char(t));
                    keyboard
					warning('bkkGenerate2:debug','Debugging... could not get positivity of %s from assumptions',char(t))
				end
			end
			
		end
		pos=collect(pos,s);
		neg=collect(neg,s);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Supporting Functon DoEvals
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function[plus minus]=DoEvals(pos,neg,params,~,phaseinterval)
	syms s ;
	minparams=cell(1,length(params));
    maxparams=cell(1,length(params));
    for i=1:length(params)
		minparams{i}=sym(['min' char(params(i))]);
		maxparams{i}=sym(['max' char(params(i))]);
    end
    
    % add the gain to the params for this last construction: % WORKING?
    params=[params 'K' ];
    minparams=[minparams 'minK'];
    maxparams=[maxparams 'maxK'];
    
    % now for phase:
    if isnumeric(phaseinterval)
        params=[params 'sin(m)' 'cos(m)'];
        
        % This is experimental at the moment. This needs further testing
        % and more robustness. I know there are some cases where this will
        % not work. For instance, in a MIMO system, one might get
        % exp(2*i*theta) as a term in the polynomial. The current
        % implementaton here will not correctly handle such a case. 
        %
        % Also, fundamentally, this provides very conservative results
        % because of the overbounding of the vertex polynomials produced
        % from the algorithm. 
        
        switch phaseinterval
            case 1
                minparams=[minparams 'sin(maxm)' 'cos(minm)'];
                maxparams=[maxparams 'sin(minm)' 'cos(maxm)'];
            
            case 2    
                minparams=[minparams 'sin(maxm)' 'cos(minm)'];
                maxparams=[maxparams 'sin(minm)' 'cos(maxm)'];
                                
            case 3
                minparams=[minparams 'sin(minm)' 'cos(maxm)'];
                maxparams=[maxparams 'sin(maxm)' 'cos(minm)'];
                
            case 4
                minparams=[minparams 'sin(maxm)' 'cos(minm)'];
                maxparams=[maxparams 'sin(minm)' 'cos(maxm)'];
            otherwise
                error('Unknown phase interval');
        end
    end
    
    %keyboard
	plus=sym(0);
	minus=sym(0);
	
	
    
    
    %get max degree of pos and neg
	degree=max(double([feval(symengine,'degree',pos,s),feval(symengine,'degree',neg,s)]) );
	for k=degree:-1:0
		this=feval(symengine,'coeff',pos,s,k); %this is pos
		that=feval(symengine,'coeff',neg,s,k); %that is neg
        
       
        % switches based on orthant
        switch double(mod(k,4))
            case 0
                plus=plus +   (subs(this,params,maxparams) + subs(that,params,minparams))*s^k;
                minus=minus + (subs(this,params,minparams) + subs(that,params,maxparams))*s^k;
            case 1
                plus=plus +   ( subs(this,params,maxparams) + subs(that,params,minparams))*s^k;
                minus=minus + ( subs(this,params,minparams) + subs(that,params,maxparams))*s^k;
            case 2
                plus=plus +   ( subs(this,params,minparams) + subs(that,params,maxparams))*s^k;
                minus=minus + ( subs(this,params,maxparams) + subs(that,params,minparams))*s^k;
            case 3
                plus=plus +   ( subs(this,params,minparams) + subs(that,params,maxparams))*s^k;
                minus=minus + ( subs(this,params,maxparams) + subs(that,params,minparams))*s^k;
            otherwise
                disp('Error in switch statement')
        end
        
        %keyboard
        
	end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Contact us