pwa2hys.m

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

function [] = pwa2hys;
%===============================================================================
%
% Title:        pwa2hys.m                                             	
%                                                                       
% Project:      Hybrid systems   
%                                                             
% Purpose:      converts pwa approximation into hysdel code
%               
% Author:       Tobias Geyer
%                                                                       
% History:      date         subject                                       
%               2001.08.02   Initial Version   
% 				2001.12.10	 switch to disable generation of bounds
%							 (so far only implemented for DaMust)
%               2002.02.05   AdDa_2reg updated
%               2002.02.12	 0*xi and 1* are ommitted
%               2002.03.01	 handles now also the case with 1 region
%  
% Usage:        The pwa approximation 
%                 y{i} = f{i}*u + g{i}  if  Hi{i}*u <= Ki{i}  for i=1..r
%               is rewritten as HYSDEL code in a way explained in the respective
%               subfunction.
%
% Remark:       Only the hyperplanes which are seperating regions are kept for HYSDEL
%               I.e. the hyperplanes which deliminate the domain of the functions are not kept
%
% Requires:     con_list6 (from Mato)
%
% Input:        see the section 'Input variables'
%
% Output:       file containing the hysdel code
%
% Contact:      Tobias Geyer
%               Automatic Control Laboratory
%               ETH Zentrum, 
%               Zurich, Switzerland
%
%               [e-mail]
%
%               Comments and bug reports are highly appreciated
%
%===============================================================================


% ------------------- Input variables (adapt this part to your needs) ----------

% the names of the input and output variables
% eg. for 2 input variables x1, x2 and 1 output variable y define the 
% following: u{1}='x1', u{2}='x2'
%            y='y'
u=[]; u{1}='U4d'; u{2}='U4q'; %u{3}='Pl_z';
y=[]; y='U4m2'; 

% number of first delta-variable to be created
first_d = 1;	

% number of first z-variable to be created
first_z = 1;

% the filename of the .mat-file in which the pwa approximation is stored
in_fname = ['/home/geyer/hybrid/ABB/approx/pwa_' y '.mat'];

% the filename of the .hys-file in which the pwa approx. is written
out_fname = ['/home/geyer/hybrid/ABB/mld/' y '.hys'];

% shift upper and lower bounds in AD and DA section by Eps
Eps = 1e-2;

% switch to disable the generation of the bounds of the AD and DA sections
% (they are not needed anymore in the new verions of HYSDEL...)
do_bounds = 0;		% 0=disabled, 1=enabled

% ------------------- don't change the rest of the file ------------------------


% ------------------- interface to Gianni's pwa approximation ------------------

% run Gianni's code
eval(['load ' in_fname]);	% load Gianni's pwa approximation
% this gives: estimated_regions: est{i,1} <= est{i,2} and
%             theta_est

%estimated_regions = estimated_regions_sim;

regions = size(estimated_regions, 1);
Hi = []; 
Ki = [];
for i = 1:regions
	Hi{i} = estimated_regions{i,1};
	Ki{i} = estimated_regions{i,2};
end;
clear estimated_regions;
% this gives: Hi{i}*x <= Ki{i}

Po = theta_est;
clear theta_est



% ------------------- test integrity of input data -----------------------------
nu = size(u, 2);

% is the number of input variables correct?
nu_nesc = size( Hi{1}, 2 );
if nu ~= nu_nesc
	error(['Number of specified input variables is wrong. Expected '...
	num2str(nu_nesc) ' input variables.']);
end;



% --------------------- preliminary stuff --------------------------------------

% run Mato's code and get neighbouring edges of hyberplanes
if regions >= 2
	[BC, Hn, Kn] = con_list6(Hi, Ki);
	% this gives e.g.: 
	% BC{1} = [2 3 0 0 0]'
	% for region #1: 1st hyperplane is adjacent to region #2
	%                2nd hypberplane is adjacent to region #3
	%                3rd, 4th and 5th hypperplanes have no neighbours
else
	BC = []; Hn = []; Kn = [];
end;

% prepare the fields
AD = [];
DA = [];
MUST = [];



% -------------------- decide which conversion function to use -----------------

if regions > 2
		
	% old approach with one delta-variable for every hyperplane
    %[AD, DA, num_d] = AdDa(Hi, Ki, Po, BC, regions, y, u, nu, first_d, first_z, Eps);
    %disp('used AD - DA approach (one delta-variable for every hyperplane)')
	
	% new approach with one delta-variable for every region (according to BeMo99)
	[DA, MUST, num_d] = DaMust(Hi, Ki, Po, BC, regions, y, u, nu, first_d, first_z, Eps, do_bounds);
    disp('used DA - MUST approach (one delta-variable for every region)')	
    
elseif regions == 2
   
   	[AD, DA, num_d] = AdDa_2reg(Hi, Ki, Po, BC, regions, y, u, nu, Eps, do_bounds);	
	disp('used AD - DA approach for 2 regions (one delta-variable)')
   
else
	
	[MUST, num_d] = one_reg(Po{1}, y, u, nu);
    disp('used MUST approach for 1 region (no delta-variable)')	
	
end;



% ----------------- store the AUX, DA, CONTINUOUS and MUST sections -------------
tab = '    ';

% open file to write HYSDEL code
fid = fopen(out_fname, 'wt');

if regions > 1

	% write AUX section
	fprintf(fid, [tab 'AUX {\n']);

	fprintf(fid, [tab tab 'BOOL ']);
	if num_d > 1
		for i = 1:num_d-1,fprintf(fid, [y '_d' num2str(first_d+i-1) ', ']); end;
		fprintf(fid, [y '_d' num2str(first_d+num_d-1) ';\n']);
	else
		fprintf(fid, [y, '_d;\n']);
	end;

	fprintf(fid, [tab tab 'REAL ']);
	if num_d == 1
		fprintf(fid, [y '_z;\n']);
	else
		fprintf(fid, [y '_z, ']);
		for i = 1:num_d-1, fprintf(fid, [y '_z' num2str(first_z+i-1) ', ']); end;
		fprintf(fid, [y '_z' num2str(first_z+num_d-1) ';\n']);
	end;

	fprintf(fid, [tab '}\n\n']);

	% write AD section
	AD_length = size(AD, 2);
	if AD_length > 0
		fprintf(fid, [tab 'AD {\n']);
		for i = 1:AD_length, fprintf(fid, [tab tab AD{i} '\n']); end;
		fprintf(fid, [tab '}\n\n']);
	end;

	% write DA section
	fprintf(fid, [tab 'DA {\n']);
	for i = 1:regions, fprintf(fid, [tab tab DA{i} '\n']); end;
	fprintf(fid, [tab '}\n\n']);

	% write CONTINUOUS section
	%fprintf(fid, [tab 'CONTINUOUS {\n']);
	%fprintf(fid, [tab tab y ' = ']);
	%if num_d == 1
	%	fprintf(fid, [y '_z;\n']);
	%else
	%	for i = 1:num_d-1, fprintf(fid, [y '_z' num2str(first_z+i-1) '+']); end;
	%	fprintf(fid, [y '_z' num2str(first_z+num_d-1) ';\n']);
	%end;
	%fprintf(fid, [tab '}\n\n']);
	
end;

% write MUST section
MUST_length = size(MUST, 2);
if MUST_length > 0
	fprintf(fid, [tab 'MUST {\n']);
	for i = 1:MUST_length, fprintf(fid, [tab tab  MUST{i} '\n']); end;
	fprintf(fid, [tab '}\n']);
end;

% close file
fclose(fid);

%
disp(['pwa approximation over ' num2str(regions) ' regions.'])
disp(['wrote ' out_fname ' successfully.'])
disp(' ')







% ==============================================================================
% ==================== subfunctions ============================================
% ==============================================================================

% -------------------- DaMust --------------------------------------------------
function [DA, MUST, num_d] = DaMust(Hi, Ki, Po, BC, regions, y, u, nu, first_d, first_z, Eps, do_bounds);
%
% regions > 2
% the conversion is done according to BeMo99 p. 411,412
%
% i.e.  AUX {
%           BOOL Y_d1, ... Y_dr;	
%           REAL Y_z, Y_z1, ... Y_zr;
%       }
%       DA {
%           Y_z1 = {IF Y_d1 THEN f{1}*un + g{1} [M1, m1, 0]};	  
%           ...
%           Y_zr = {IF Y_dr THEN f{r}*un + g{r} [Mr, mr, 0]};
%       }
%       MUST {
%           Hi{1}*u - Ki{1} - N1 + N1*d1 <= 0;
%           ...
%           Hi{r}*u - Ki{r} - Nr + Nr*dr <= 0;
%           Y_d1 + ... + Y_dr = 1;
%           Y_z = Y_z1 + ... + Y_zr;
%       }
%
% where r = number of regions

% create DA entries
for r = 1:regions
	hi = Hi{r}; ki = Ki{r};
	po = Po{r};
	
	% add ... = {IF .. THEN 
	da = [y '_z' num2str(first_z+r-1) ' = {IF ' y '_d' num2str(first_d+r-1) ' 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))];
		
	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;
		
	% finish line
	da = [da '};'];
	
	% store new DA line
	DA{r} = da;	
end;


% create MUST section
m = 1;	% line counter of MUST section
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
			% we have found a hyperplane seperating two neighbouring regions
			% we will rewrite this hyperplane into the MUST section
			
			MUST{m} = '';
			
			% ex.: '-6.3511*I2d -2.6294*I2q '
			first = 1;
			for i = 1:nu
				if hi(c,i) ~= 0
					if (hi(c,i) >= 0) & ~first, MUST{m} = [MUST{m} '+']; end;
					if hi(c,i) == (-1), 
						MUST{m} = [MUST{m} '-'];
					elseif hi(c,i) ~= 1, 
						MUST{m} = [MUST{m} num2str(hi(c,i)) '*']; 
					end;
					MUST{m} = [MUST{m} u{i} ' '];
					first = 0;
				end;
			end;

			% ex.: '+5.0265'
			if -ki(c) >= 0, MUST{m} = [MUST{m} '+']; end;
			MUST{m} = [MUST{m} num2str(-ki(c)) ' '];
								
			% get upper 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_max = hp_max + Eps;

			% ex.: ' -100'
			if -hp_max >= 0, MUST{m} = [MUST{m} '+']; end;
			MUST{m} = [MUST{m} num2str(-hp_max) ' '];
			
			% ex.: ' +100*(REAL d2) <= 0;'
			if hp_max >= 0, MUST{m} = [MUST{m} '+']; end;
			MUST{m} = [MUST{m} num2str(hp_max) '*(REAL ' y '_d' num2str(first_d+r-1) ') <= 0;'];
			
			% increase the line counter of the MUST section
			m = m+1;
		
		end;
	end;
		
end;

% create equality constraints in MUST section
% first: d1+d2+...dn - 1 <= 0;
MUST{m} = '';
for i = first_d:first_d+regions-2
	MUST{m} = [MUST{m} '(REAL ' y '_d' num2str(i) ') +'];
end;
MUST{m} = [MUST{m} '(REAL ' y '_d' num2str(i+1) ') -1 <= 0;'];
m = m+1;

% then: -d1-d2-...-dn + 1 <= 0;