pwa2hys.m

由matlab开发的hybrid系统的描述语言

MUST{m} = '';
for i = first_d:first_d+regions-1
	MUST{m} = [MUST{m} '-(REAL ' y '_d' num2str(i) ') '];
end;
MUST{m} = [MUST{m} '+1 <= 0;'];
m = m+1;

% create additional equality constraints in MUST section
% first  y_z -(y_z1+y_z2+...y_zm) <= 0;
MUST{m} = [' ' y '_z -(' y '_z1'];
for i = first_d+1:first_d+regions-1
	MUST{m} = [MUST{m} '+' y '_z' num2str(i)];
end;
MUST{m} = [MUST{m} ') <= 0;'];
m = m+1;

% then: -y_z + y_z1+yz_z2+...+y_zm  <= 0;
MUST{m} = ['-' y '_z  '];
for i = first_d:first_d+regions-1
	MUST{m} = [MUST{m} '+' y '_z' num2str(i)];
end;
MUST{m} = [MUST{m} '  <= 0;'];
m = m+1;

% set number of delta variables
num_d = regions;


return;




% -------------------- AdDa_2reg -----------------------------------------------
function [AD, DA, num_d] = AdDa_2reg(Hi, Ki, Po, BC, regions, y, u, nu, Eps, do_bounds);
% regions = 2
% 
% i.e.  AUX {
%           BOOL Y_d;	
%           REAL Y_z;
%       }
%       AD {
%           Y_d = Hi*u - Ki <= 0 [M, m, eps];
%       }
%       DA {
%           Y_z = {IF d1 THEN f{1}*u + g{1} [N2, n2, 0]
%                        ELSE f{2}*u + g{2} [N2, n2, 0]};
%       }
%       CONTINUOUS {
%           Y = Y_z;	
%       }

options = optimset('DISPLAY', 'off');


% AD section -----------------------

% Look at the first region
bc_curr = BC{1};
hi = Hi{1};
ki = Ki{1};

% Find the hyperplane separating region #1 and region #2
c = 1;
while (bc_curr(c) == 0) & (c <= length(bc_curr))
	c = c+1;
end;
% ==> hp = hi(c,:) - ki(c) is the hyperplane we have been looking for

% ex.: 'y_d = '
ad = [y '_d = '];

% ex.: '-6.3511*I2d -2.6294*I2q '
first = 1;
for i = 1:nu
	if hi(c,i) ~= 0
		if (hi(c,i) >= 0) & ~first, ad = [ad '+']; end;
		if hi(c,i) == (-1),
			ad = [ad '-']; 
		elseif hi(c,i) ~= 1
			ad = [ad num2str(hi(c,i)) '*']; 
		end;
		ad = [ad u{i} ' '];
		first = 0;
	end;
end;

% ex.: '+5.0265 <= 0'
if -ki(c) >= 0, ad = [ad '+']; end;
ad = [ad num2str(-ki(c)) ' <= 0'];

if do_bounds
	% get upper and lower bound 
	% i.e. evaluate hyperplane hp = hi(c,:) - ki(c) over all regions
	[hp_max, hp_min] = bounds(hi(c,:), -ki(c), Hi, Ki, regions);
	hp_min = hp_min - Eps;
	hp_max = hp_max + Eps;

	% ex.: ' [100, -100, eps];'
	ad = [ad ' [' num2str(hp_max) ', ' num2str(hp_min) ', eps];'];
else
	ad = [ad ';'];
end;

% store AD{1}
AD{1} = ad;


% DA section -------------------------

% prepare the DA fields:
% DA{1} = 'Yz = {IF d1 THEN '
% DA{2} = '            ELSE '
DA{1} = [y '_z = {IF ' y '_d '];
DA{2} = '';
for i = 1:length(DA{1})		% fill DA{2} with spaces
	DA{2} = [DA{2} ' '];
end;
DA{1} = [DA{1} 'THEN '];
DA{2} = [DA{2} 'ELSE '];

% add polynom and bounds to DA entries
for r = 1:2
	da = DA{r};
	hi = Hi{r}; ki = Ki{r};
	po = Po{r};
	
	% add polynom
	for i = 1:nu
		if (po(i) >= 0) & (i > 1), da = [da '+']; end;
		da = [da num2str(po(i)) '*' u{i} ' '];
	end;
	if po(i+1) >= 0, da = [da '+']; end;
	da = [da num2str(po(i+1))];
	
	if do_bounds
		% get upper and lower bound 
		% i.e. evaluate pwa approximation po over all regions
		[up_bound, lo_bound] = bounds(po(1:nu), po(nu+1), Hi, Ki, regions);
		up_bound = up_bound + Eps;
		lo_bound = lo_bound - Eps;

		% add bounds
		da = [da ' [' num2str(up_bound) ', '];
		da = [da num2str(lo_bound) ', 0]'];
	end;
	
	% add '};' if necessary
	if r == 2, da = [da '};']; end;

	% store new DA line
	DA{r} = da;	
end;

% set number of delta variables
num_d = 1;

return




% -------------------- one_reg -------------------------------------------------
function [MUST, num_d] = one_reg(po, y, u, nu);
%
% regions = 1
%
%       MUST {
%           Y = f*un + g;
%       }
%

% first:  y -(po{1}*u_1 +..+ po{nu}*u_nu + po{nu+1}) <= 0;
MUST{1} = [' ' y ' -('];
for i = 1:nu
	if (po(i) >= 0) & (i > 1), MUST{1} = [MUST{1} '+']; end;
	MUST{1} = [MUST{1} num2str(po(i)) '*' u{i} ' '];
end;
if po(i+1) >= 0, MUST{1} = [MUST{1} '+']; end;
MUST{1} = [MUST{1} num2str(po(i+1)) ') <= 0;'];

% then: -y + po{1}*u_1 +..+ po{nu}*u_nu + po{nu+1}  <= 0;
MUST{2} = ['-' y '  '];
for i = 1:nu
	if (po(i) >= 0) & (i >= 1), MUST{2} = [MUST{2} '+']; end;
	MUST{2} = [MUST{2} num2str(po(i)) '*' u{i} ' '];
end;
if po(i+1) >= 0, MUST{2} = [MUST{2} '+']; end;
MUST{2} = [MUST{2} num2str(po(i+1)) '  <= 0;'];

% set number of delta variables
num_d = 1;


return;




% -------------------- AdDa ----------------------------------------------------
function [AD, DA, num_d] = AdDa(Hi, Ki, Po, BC, regions, y, u, nu, first_d, fist_z, Eps);
% regions > 2
% introduce a delta-variable for every hyperplane seperating two adjacent regions
% i.e.  AUX {
%           BOOL d1, ... dm;	
%           REAL Yz1, ... Yzr;
%       }
%       AD {
%           d1 = Hi{i}(j)*u1 + ... + Hi{i}(j)*un - Ki{j} <= 0 [?, -?, eps];
%           ...
%          dm = ...
%       }
%       DA {
%          Yz1 = {IF d? & d? THEN f{i}(1)*u1 + ... +f{i}(n)*un + g{i} [?, -?, 0]};	  
%          ...
%          Yzr = ...
%       }
%       CONTINUOUS {
%          Y = Yz1 + ... + Yzr;	
%       }

% prepare the DA field
for r = 1:regions
	DA{r} = [y 'z' num2str(first_z+r-1) ' = {IF '];
end;

options = optimset('DISPLAY', 'off');
d = first_d;

for r = 1:regions
	bc_curr = BC{r};
	hi = Hi{r};
	ki = Ki{r};

	for c = 1:length(bc_curr)
		if ~bc_curr(c) == 0
			s = bc_curr(c);			% r = # of current region
									% s = # of neighbouring region
			bc_neigh = BC{s};
			
			for n = 1:length(bc_neigh)
				if bc_neigh(n) == r
					% we have found the corresponding hyperplanes
					
					% create new AD entry (the d-th one):
					% e.g.: d1 = -6.3511*I2d -2.6294*I2q +5.0265 <= 0 [100, -100, eps];
		
					% ex.: 'd1 = '
					AD{d} = ['d' num2str(d) ' = '];
					
					% ex.: '-6.3511*I2d -2.6294*I2q '
					for i = 1:nu
						if (hi(c,i) >= 0) & (i > 1), AD{d} = [AD{d} '+']; end;
						AD{d} = [AD{d} num2str(hi(c,i)) '*' u{i} ' '];
					end;
					
					% ex.: '+5.0265'
					if -ki(c) >= 0, AD{d} = [AD{d} '+']; end;
					AD{d} = [AD{d} num2str(-ki(c))];
															
					% get upper and lower bound 
					% i.e. evaluate hyperplane hp = hi(c,:) - ki(c) over all regions
					[hp_max, hp_min] = bounds(hi(c,:), -ki(c), Hi, Ki, regions);
					hp_min = hp_min - Eps;
					hp_max = hp_max + Eps;
						
					% ex.: ' <= 0 [100, -100, eps];'
					AD{d} = [AD{d} ' <= 0 [' num2str(hp_max) ', ' num2str(hp_min) ', eps];'];
										
					% update DA entries:
					DA{r} = [DA{r}  'd' num2str(d) ' & '];
					DA{s} = [DA{s} '~d' num2str(d) ' & '];
										
					% we set the respective entries of bc to zero:
					bc_curr(c) = 0;
					bc_neigh(n) = 0;
					
					% increase the counter of the delta variables
					d = d+1;
				end;
			end;
			
			BC{s} = bc_neigh;			
		end;
	end;
			
	BC{r} = bc_curr;
end;

% add polynom and bounds to DA entries
for r = 1:regions
	da = DA{r};
	hi = Hi{r}; ki = Ki{r};
	po = Po{r};
	
	% delete last '&'
	da = da(1:length(da)-3);
	
	% add 'THEN'
	da = [da ' THEN '];
	
	% add polynom
	for i = 1:nu
		if (po(i) >= 0) & (i > 1), da = [da '+']; end;
		da = [da num2str(po(i)) '*' u{i} ' '];
	end;
	if po(i+1) >= 0, da = [da '+']; end;
	da = [da num2str(po(i+1))];
	
	% get upper and lower bound 
	% i.e. evaluate pwa approximation po over all regions
	[up_bound, lo_bound] = bounds(po(1:nu), po(nu+1), Hi, Ki, regions);
	up_bound = up_bound + Eps;
	lo_bound = lo_bound - Eps;
	
	% add bounds
	da = [da ' [' num2str(up_bound) ', '];
	da = [da num2str(lo_bound) ', 0]};'];

	% store new DA line
	DA{r} = da;	
end;

% set number of delta variables
num_d = d-1;

return;





% ==============================================================================
% ==================== subsubfunctions =========================================
% ==============================================================================

function [up_bound, lo_bound] = bounds(f, g, Hi, Ki, regions);
% get upper and lower bounds
% i.e. maximize and minimize f*x+g over all regions.
% these regions are defined by Hi{i}*x <= Ki{i}, i=1..regions
	
options = optimset('DISPLAY', 'off');

for t = 1:regions
	% min
	[x, fval, exitflag] = linprog(f, Hi{t}, Ki{t}, [],[],[],[],[], options);
	if exitflag <= 0, warning('unable to determine lower bound'); end;
	v = f*x + g;

	% max
	[x, fval, exitflag] = linprog(-f, Hi{t}, Ki{t}, [],[],[],[],[], options);
	if exitflag <= 0, warning('unable to determine upper bound'); end;
	w = f*x + g;

	% update minimum and maximum
	if t == 1
		lo_bound = v;
		up_bound = w;
	else
		lo_bound = min( lo_bound, v );
		up_bound = max( up_bound, w );
	end;
end;

return;