% Copyright c 1998-2002 The Board of Trustees of the University of Illinois
% 		  All rights reserved.
% Developed by:	Large Scale Systems Research Laboratory
%               Professor Richard Braatz, Director
%               Department of Chemical Engineering
%		University of Illinois
%		http://brahms.scs.uiuc.edu
% 
% Permission hereby granted, free of charge, to any person obtaining a copy
% of this software and associated documentation files (the "Software"), to
% deal with the Software without restriction, including without limitation
% the rights to use, copy, modify, merge, publish, distribute, sublicense,
% and/or sell copies of the Software, and to permit persons to whom the 
% Software is furnished to do so, subject to the following conditions:
% 		1. Redistributions of source code must retain the above copyright
%		   notice, this list of conditions and the following disclaimers.
%		2. Redistributions in binary form must reproduce the above 
%		   copyright notice, this list of conditions and the following 
%		   disclaimers in the documentation and/or other materials 
%		   provided with the distribution.
%		3. Neither the names of Large Scale Research Systems Laboratory,
%		   University of Illinois, nor the names of its contributors may
%		   be used to endorse or promote products derived from this 
%		   Software without specific prior written permission.
%
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
% OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
% FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL 
% THE CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
% OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, 
% ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
% DEALINGS IN THE SOFTWARE.
%
% SR Successive Multidimensional Realization 
%
% 	[Lr,bsr] = sr(L,bs)
%
% 	This produces a reduced uncertain systems for uncertain systems
%	written in the form of LFT (Linear Fractional Transformation) 
% 	with repeated scalar uncertainties.
%	
%	[Lr,bsr] = sr(L,bs,tol) specifies the tolerance used to
%	eliminate uncontrollable and unobservable subspaces. The
%	tolerance is exactly that in minreal. 
%
%	Inputs: L - system matrix, CONSTANT
%		bs - block structure of L, which includes the dimension
%		     of the output as the last element in bs
%
%	Outputs: Lr - system matrix for reduced system, CONSTANT
%		 bsr - block structure of Lr, which includes the dimension
%		     of the output as the last element in bsr
%
%	Algorithm:  The Successive Multidimensional Realization algorithm
%	is based on repeated calls to a one-dimensional realization algorithm.
%	The implementation below is based on minreal because that program is
%	a part of the Control System Toolbox, which is available to most
%	users of MATLAB.  Users with access to the mu-analysis and synthesis 
%	toolbox should replace the use of minreal.m with sysbal.m to provide
%	enhanced numerical robustness (the general principles behind this
%	are described in the 1D case by Bruce Moore in his 1981 IEEE TAC paper).
%
%
%	Users of this code should cite the following journal papers:
%
%	E.L. Russell, C.P.H. Power, and R.D. Braatz. Multidimensional
%	realizations of large scale uncertain systems for multivariable
%	stability margin computation, Int. J. of Robust and Nonlinear
%	Control, 7:113-125, 1997.
%
%	R. Gunawan, E.L. Russell, and R.D. Braatz. Comparison of theoretical
%	and computational characteristics of dimensionality reduction methods
%	for large scale uncertain system, Journal of Process Control, 
%	11:543-552, 2001.
%
% Authors       : Evan L. Russell, Rudiyanto Gunawan, and Richard D. Braatz
% Last Modified : January 2 20, 2002
% Email         : braatz@uiuc.edu
% WWW	        : http://brahms.scs.uiuc.edu/lssrl/software/

function [Lr,bsr]=sr(L,bs,tol);

r = length(bs);
bsr = bs;
Lr = L;
red = 0;

for i=1:r-1
	[nr,nc]=size(Lr);
        cumbs=[0 cumsum(bsr)];
	ni=cumbs(i)+1:cumbs(i+1);
        ni1=1:cumbs(i); nri2=cumbs(i+1)+1:nr; nci2=cumbs(i+1)+1:nc;
        A=Lr(ni,ni);, B=Lr(ni, [ni1 nci2]);
	C=Lr([ni1 nri2],ni);, D=Lr([ni1 nri2],[ni1 nci2]);

	if (exist('tol')==0),
		[Am,Bm,Cm,Dm]=minreal(A,B,C,D);
	else
		[Am,Bm,Cm,Dm]=minreal(A,B,C,D,tol);
	end

	Lr=[Am Bm;Cm Dm];
	ri=size(Am,1); redi=bsr(i)-ri; red=red+redi; bsr(i)=ri;  
        nr=[ri+ni1 1:ri nri2-redi]; nc=[ri+ni1 1:ri nci2-redi];
        Lr=Lr(nr,nc);
end