#! /bin/sh # This is a shell archive. Remove anything before this line, then unpack # it by saving it into a file and typing "sh file". To overwrite existing # files, type "sh file -c". You can also feed this as standard input via # unshar, or by typing "sh 'READ.ME' <<'END_OF_FILE' X *************************************************************************** X * All the software contained in this library is protected by copyright. * X * Permission to use, copy, modify, and distribute this software for any * X * purpose without fee is hereby granted, provided that this entire notice * X * is included in all copies of any software which is or includes a copy * X * or modification of this software and in all copies of the supporting * X * documentation for such software. * X *************************************************************************** X * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED * X * WARRANTY. IN NO EVENT, NEITHER THE AUTHORS, NOR THE PUBLISHER, NOR ANY * X * MEMBER OF THE EDITORIAL BOARD OF THE JOURNAL "NUMERICAL ALGORITHMS", * X * NOR ITS EDITOR-IN-CHIEF, BE LIABLE FOR ANY ERROR IN THE SOFTWARE, ANY * X * MISUSE OF IT OR ANY DAMAGE ARISING OUT OF ITS USE. THE ENTIRE RISK OF * X * USING THE SOFTWARE LIES WITH THE PARTY DOING SO. * X *************************************************************************** X * ANY USE OF THE SOFTWARE CONSTITUTES ACCEPTANCE OF THE TERMS OF THE * X * ABOVE STATEMENT. * X *************************************************************************** X X AUTHORS: X X W. KRAJEWSKI X POLISH ACADEMY OF SCIENCE - POLAND X X A. LEPSCHY, M. REDIVO-ZAGLIA, U. VIARO X UNIVERSITY OF PADOVA - ITALY X X X REFERENCE: X X - A PROGRAM FOR SOLVING THE L_2 REDUCED-ORDER MODEL PROBLEM WITH X FIXED DENOMINATOR DEGREE X NUMERICAL ALGORITHMS, 9(1995), PP. 355-377 X X SOFTWARE REVISION DATE: X X MAY 24, 1995 X X SOFTWARE LANGUAGE: X X MATLAB X END_OF_FILE if test 1820 -ne `wc -c <'READ.ME'`; then echo shar: \"'READ.ME'\" unpacked with wrong size! fi # end of 'READ.ME' fi if test -f 'energy.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'energy.m'\" else echo shar: Extracting \"'energy.m'\" \(1017 characters\) sed "s/^X//" >'energy.m' <<'END_OF_FILE' Xfunction err2 = energy(ner,der) X% X% Description: X% This function computes the squared norm err2 of X% E(s) = ner(s)/der(s) (error between the original and the X% reduced transfer functions) according to the Astrom method X% (see "Introduction to Stochastic Control"). X% X% Usage: X% err2 = energy(ner,der); X% X% Input parameters: X% ner - numerator polynomial of the error function (descending X% powers of s) X% der - denominator polynomial of the error function (descending X% powers of s) X% X% Output parameter: X% err2 - squared error norm (index value) X% X Xa = der; Xb = ner; Xn = length(a) - 1; Xif length(b) < n X b = [zeros(1,n - length(b)),b]; Xend Xv = 0; Xif a(1) <= 0 X err2 = 0; X return Xend Xfor i = 1:n X if a(i+1) <= 0 X err2 = 0; X return X end X alfa = a(i)/a(i+1); X beta = b(i)/a(i+1); X v = v + (beta/alfa)*beta; X i1 = i + 2; X if i1 - n <= 0 X for j = i1:2:n X a(j) = a(j) - alfa*a(j+1); X b(j) = b(j) - beta*a(j+1); X end X end Xend Xerr2 = v/2; END_OF_FILE if test 1017 -ne `wc -c <'energy.m'`; then echo shar: \"'energy.m'\" unpacked with wrong size! fi # end of 'energy.m' fi if test -f 'l2mimo.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'l2mimo.m'\" else echo shar: Extracting \"'l2mimo.m'\" \(7978 characters\) sed "s/^X//" >'l2mimo.m' <<'END_OF_FILE' Xfunction [nr,dr,jr,kr,np,dp,jp,kp] = l2mimo(num,den,mi,mo,epsa,kmax,kdisp) X% X% Description: X% l2mimo performs L_2 model reduction starting from a transfer matrix X% F(s) = N(s)/den(s) with mi inputs and mo outputs via the iterative X% interpolation method described in the paper of Krajewski et al. X% N(s) is the matrix of numerators and den(s) is the common denominator. X% num is the vector formed from the coefficients of all elements of N(s) X% taken in lexicographic order. X% If deg den(s) = ni, then the elements of N are polynomials of degree at X% most ni-1. The subvectors of num corresponding to these polynomials must X% consist of ni elements some of which may be zero (e.g. when the degree of X% a given numerator is less than ni-1). In particular, it may happen that X% the whole subvector consists of zeros only. Thus, num must be an mi*mo*ni X% vector. Coefficients of polynomials are arranged according to ascending X% powers of s. X% X% The procedure l2mimo uses Matlab fliplr and flipud commands. X% In some versions of Matlab these commands are not available X% and must be replaced by flipy and flipx. X% X% The procedure l2mimo uses the following functions: X% energy to evaluate the index X% minr to input the data for the reduced model X% mour to display the final solution X% X% Usage: X% [nr,dr,jr,kr,np,dp,jp,kp] = l2mimo(num,den,mi,mo,epsa,kmax,kdisp); X% X% Input parameters: X% num - vector of numerators coefficients stored according to the X% above Description. X% If the element (i,j) of the numerator matrix N(s) equals X% n(i,j)(s) = n(i,j,0)+n(i,j,1)*s +...+ n(i,j,ni-1)*s^(ni-1) X% for i = 1,...,mo, j = 1,...,mi, then the vector num is as follows X% num = [n(1,1,0) ... n(1,1,ni-1) n(1,2,0) ... n(1,2,ni-1) ... X% n(mo,mi,0) ... n(mo,mi,ni-1)] X% All coefficients, including zero coefficients, must be specified. X% den - original model denominator (monic) X% If the denominator takes the form X% den(s) = d(0) + d(1)*s + ... + d(ni-1)*s^(ni-1) + s^(ni), X% vector den is as follows X% den = [d(0) d(1) ... d(ni-1) d(ni)], with d(ni)=1. X% mi - number of inputs X% mo - number of outputs X% Optional: X% epsa - accuracy of the approximation; if not specified epsa = 0.001 X% kmax - maximum number of iterations; if not specified kmax = 50 X% kdisp - if given, every kdisp iterations partial results (denominator X% coefficients) are displayed X% X% Output parameters: X% nr,dr,jr,kr - reduced order numerators, denominator, value of the X% minimized index and number of iterations when the stopping X% criterion is satisfied. X% np,dp,jp,kp - reduced order numerators, denominator, value of the X% minimized index and number of iterations corresponding to the X% least value of the index within the kr or kmax iterations (jp X% may be different from the value jr of the index corresponding to X% the satisfaction of the chosen stopping criterion). X% X X% X% ni - order of the original model X% Xni = length(den) - 1; X% X% Specification of the reduced model order and first guesses X% X% dg0 - first guess of denominator. X% ng0 - first guess of numerators. X% X[ng0,dg0] = minr(mi,mo,ni); X% X% r - order of the reduced model X% Xr = length(dg0) - 1; X% X% m - total number of elements of the numerator vector (equal to mi*mo*ni) X% Xm = length(num); X% X% ito - number of the polynomial entries n(i,j) of the numerator X% matrix N(s) X% Xito = mi*mo; X% X% Default values X% Xif nargin < 5 X epsa = 0.001; X kmax = 50; X kdisp = kmax + 1; Xelseif nargin == 5 X kmax = 50; X kdisp = kmax + 1; Xelseif nargin == 6 X kdisp = kmax + 1; Xend X% X% Matrices df1, df2 X% Xvr = [den(1) zeros(1,r-1)]; Xvc = [den';zeros(r-1,1)]; Xmatp = toeplitz(vc,vr); Xdf1 = -matp(1:r,:); Xdf2 = -matp(r+1:ni+r,:); Xdf1 = eye(r)/df1; Xdf2 = df2*df1; X% X% Initializations X% Xdg = (dg0(1:r))'; Xng = ng0; Xk = 0; Xsw = 1; Xdelta = inf; Xdgprev = dg; Xksw = kdisp; Xformat long; Xkp = 0; Xjp = 0; Xjg = 1; Xjf = 1; Xdp = dg0; Xnp = ng0; X% X% jp is evaluated using the Routh-like algorithm due to Astrom X% (see function energy) X% Xfor i = 1:ito X jp = jp + energy(conv(fliplr(num(jf:jf+ni-1)),fliplr(dg0)) - ... X conv(fliplr(ng0(jg:jg+r-1)),fliplr(den)), ... X conv(fliplr(den),fliplr(dg0))); X jg = jg + r; X jf = jf + ni; Xend Xif jp == 0 X jp = inf; Xend X% X% Iterations X% Xwhile delta > epsa*min(abs(dgprev)) X if k > kmax X sw = 0; X break X end X% X% Partial results are displayed every kdisp iterations X% X if k == ksw X disp(' ') X fprintf(' At iteration %3.0f\n',k) X disp(' Partial solution (in ascending powers of s): ') X ddg = [dg;1]'; X for i = 1:r+1 X fprintf(' %18.14g',ddg(i)) X end X ksw = ksw + kdisp; X end X dgprev = dg; X ngprev = ng; X k = k + 1; X vr = [-dgprev(1) zeros(1,ni-1)]; X vc = [dgprev;1;zeros(ni-1,1)]; X ii = 1; X for i = 1:r+1 X ii = -ii; X vc(i) = ii*vc(i); X end X matp = toeplitz(vc,vr); X dg1 = matp(1:r,:); X dg2 = matp(r+1:ni+r,:); X d1 = matp(2:r+1,2:ni); X d2 = matp(r+2:ni+r,2:ni); X d1 = d1/d2; X dg2 = eye(ni)/(dg2 - df2*dg1); X aa = zeros(r,r); X cc = zeros(r,1); X jf = 1; X jg = 1; X for i = 1:ito X vr = [ngprev(jg) zeros(1,ni-1)]; X vc = [(ngprev(jg:jg+r-1))';zeros(ni-1,1)]; X ii = -1; X for j = 1:r X ii = -ii; X vc(j) = ii*vc(j); X end X matp = toeplitz(vc,vr); X ng1 = matp(1:r,:); X ng2 = matp(r+1:ni+r-1,:); X ng1 = (ng1 - d1*ng2)*dg2; X vr = [num(jf) zeros(1,r-1)]; X vc = [(num(jf:jf+ni-1))';zeros(r,1)]; X matp = toeplitz(vc,vr); X nf1 = matp(1:r,:); X nf2 = matp(r+1:ni+r,:); X aa = aa + ng1*(nf2 - df2*nf1); X ci = -(num(jf:jf+ni-1))'; X cc = cc + ng1*ci; X jg = jg + r; X jf = jf + ni; X end X dg = aa\cc; X% X% Index value for the solution obtained after k iterations X% X jg = 1; X jf = 1; X jpom = 0; X for i = 1:ito X vr = [num(jf) zeros(1,r-1)]; X vc = [(num(jf:jf+ni-1))';zeros(r,1)]; X matp = toeplitz(vc,vr); X nf1 = matp(1:r,:); X nf2 = matp(r+1:ni+r,:); X ci = -(num(jf:jf+ni-1))'; X x2 = dg2*(ci - (nf2 - df2*nf1)*dg); X x1 = -df1*(dg1*x2 + nf1*dg); X for j = 1:r X ng(jg+j-1) = x1(j); X end X jpom = jpom + energy(conv(fliplr(num(jf:jf+ni-1)),[1,fliplr(dg')]) - ... X conv(fliplr(ng(jg:jg+r-1)),fliplr(den)), ... X conv(fliplr(den),[1,fliplr(dg')])); X jf = jf + ni; X jg = jg + r; X end X% X% Comparison of the current solution with the best one obtained X% within k-1 iterations. If the current solution is better, then X% it is stored as the best solution after k iterations X% X if jpom > 0 X if jpom < jp X kp = k; X jp = jpom; X% X% The best solution is saved X% X np = ng; X dp = [dg' 1]; X end X end X delta1 = max(abs(dg - dgprev)); X delta2 = max(abs(ng - ngprev)); X delta = max(delta1,delta2); Xend X% X% After a solution has been reached, the resulting polynomials coefficients X% are displayed according to ascending powers of s X% Xif sw == 1 X kr = k; X jr = jpom; X nr = ng; X dr = [dg' 1]; X disp(' ') X fprintf(' Solution found in %3.0f iterations',k) X disp(' ') X disp(' Strike any key to continue ... '), pause X mour(nr,dr,jr,kr,mi,mo); X if kp ~= kr X disp(' A solution characterized by an index value lower than') X disp(' the one corresponding to the satisfaction of the ') X disp(' stopping criterion has been found') X disp(' Strike any key to continue ... '), pause X mour(np,dp,jp,kp,mi,mo); X end Xelse X disp(' ') X fprintf(' Cannot find any solution within %3.0f iterations',kmax) X disp(' ') X disp(' The best solution with respect to the index value follows') X disp(' Strike any key to continue ... '), pause X mour(np,dp,jp,kp,mi,mo); Xend END_OF_FILE if test 7978 -ne `wc -c <'l2mimo.m'`; then echo shar: \"'l2mimo.m'\" unpacked with wrong size! fi # end of 'l2mimo.m' fi if test -f 'l2siso.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'l2siso.m'\" else echo shar: Extracting \"'l2siso.m'\" \(6809 characters\) sed "s/^X//" >'l2siso.m' <<'END_OF_FILE' Xfunction [nr,dr,jr,kr,np,dp,jp,kp] = l2siso(num,den,epsa,kmax,kdisp) X% X% Description: X% l2siso performs L_2 model reduction starting from a transfer function X% f(s) = num(s)/den(s) via the iterative interpolation method described X% in the paper by Krajewsky et al. X% num is the original numerator and den is the original denominator. X% If deg den(s) = ni, then num(s) ia a polynomial of degree at X% most ni-1. Thus num(s) consists of ni coefficients (some of them, X% however, may be zero). X% The coefficients of the polynomials are ordered according to X% ascending powers of s. X% X% The procedure l2siso uses Matlab fliplr command. X% In some versions of Matlab this command is not available X% and must be replaced by flipy. X% X% The procedure l2siso uses the following functions: X% energy to evaluate the index X% sinr to input the data for the reduced model X% sour to display the final solution X% X% Usage: X% [nr,dr,jr,kr,np,dp,jp,kp] = l2siso(num,den,epsa,kmax,kdisp); X% X% Input parameters: X% num - original numerator X% If the numerator takes the form X% num(s) = n(0) + n(1)*s + ... + n(ni-2)*s^(ni-2) + n(ni-1)s^(ni-1) X% vector num is as follows X% num = [n(0) n(1) ... n(ni-2) n(ni-1)] X% All coefficients, including zero coefficients, must be specified. X% den - original model denominator (monic) X% If the denominator takes the form X% den(s) = d(0) + d(1)*s + ... + d(ni-1)*s^(ni-1) + s^(ni) X% vector den is as follows X% den = [d(0) d(1) ... d(ni-1) d(ni)], with d(ni)=1 X% Optional: X% epsa - accuracy of the approximation; if not specified epsa = 0.001 X% kmax - maximum number of iterations; if not specified kmax = 50 X% kdisp - if given, every kdisp iterations partial results (denominator X% coefficients) are displayed X% X% Output parameters: X% nr,dr,jr,kr - reduced order numerator, denominator, value of the X% minimized index and number of iterations when the stopping X% criterion is satisfied. X% np,dp,jp,kp - reduced order numerator, denominator, value of the X% minimized index and number of iterations corresponding to the X% least value of the index within the kr or kmax iterations (jp X% may be different from the value jr of the index corresponding to X% the satisfaction of the chosen stopping criterion). X% X X% X% ni - order of the original model X% Xni = length(den) - 1; X% X% Specification of the reduced model order and first guess X% X% dg0 - first guess for the denominator. X% Xdg0 = sinr(ni); X% X% r - order of the reduced model X% Xr = length(dg0) - 1; X% X% Default values X% Xif nargin < 3 X epsa = 0.001; X kmax = 50; X kdisp = kmax + 1; Xelseif nargin == 3 X kmax = 50; X kdisp = kmax + 1; Xelseif nargin == 4 X kdisp = kmax + 1; Xend X% X% The order of the coefficients is reversed X% Xnum = fliplr(num); Xden = fliplr(den); Xdg0 = fliplr(dg0); X% X% Submatrices m11, m21, m31 X% Xv1 = zeros(1,r); Xv2 = [0;num';zeros(r-1,1)]; Xmw = toeplitz(v2,v1); Xm11 = mw(1:ni-r,:); Xm21 = mw(ni-r+1:ni,:); Xm31 = mw(ni+1:ni+r,:); X% X% Submatrices m12, m22, m32 X% Xv1(1) = 1; Xv2 = [den';zeros(r-1,1)]; Xmw = toeplitz(v2,v1); Xm12 = -mw(1:ni-r,:); Xm22 = -mw(ni-r+1:ni,:); Xm32 = -mw(ni+1:ni+r,:); X% X% Computation of the numerator coefficients corresponding to X% the optimal denominator X% Xs12 = mw(1:ni,:); Xs22 = mw(ni+1:ni+r,:); X% X% Sub-blocks of matrix H X% Xh33 = eye(r); Xh33 = h33/m32; Xh23 = -m22*h33; Xif det(m22*h33*m31-m21) == 0 X disp(' ') X disp('Matrix m22*h33*m31-m21 is singular') X disp('Sorry, this version of l2siso cannot go further') X return Xend Xh12 = (m11-m12*h33*m31)/(m22*h33*m31-m21); Xh13 = -(m12+h12*m22)*h33; X% X% Sub-blocks of matrix N (they do not vary from one iteration to X% another) X% Xn21 = eye(r); Xn21 = n21/(m21+h23*m31); Xn31 = h33*m31; X% X% Right-hand side vectors X% Xe1 = -num(1:ni-r); Xe1 = e1'; Xe2 = -num(ni-r+1:ni); Xe2 = e2'; Xe1 = e1 + h12*e2; X% X% Determination of a candidate for the global minimum X% Xv1 = [1,zeros(1,ni-1)]; Xv2 = [dg0';zeros(ni-1,1)]; Xfor i = ni+1:2:ni+r X v2(i) = -v2(i); Xend Xv1(1) = v2(1); Xmw = toeplitz(v2,v1); Xs11 = mw(1:ni,:); Xs21 = mw(ni+1:ni+r,:); Xve1 = [m11;m21]*(dg0(2:r+1))' + num'; Xve2 = m31*(dg0(2:r+1))'; Xs21 = s21/s11; Xnp = (s22 - s21*s12)\(ve2 - s21*ve1); Xnp = np'; Xdp = dg0; X% X% jp is evaluated using the Routh-like algorithm due to Astrom X% (see function energy) X% Xjpom = energy(conv(num,dp)-conv(np,den),conv(den,dp)); Xnp = fliplr(np); Xdp = fliplr(dp); Xif jpom > 0 X jp = jpom; Xelse X jp = inf; Xend X% Xkp = 0; Xak = dg0(2:r+1); Xak = ak'; Xk = 0; Xok = 1; Xdelta = inf; Xaprev = ak; Xksw = kdisp; Xformat long; X% X% Iterations X% Xwhile delta > epsa*min(abs(aprev)) X if k > kmax X ok = 0; X break X end X% X% Partial results are displayed every kdisp iterations X% X if k == ksw X disp(' ') X fprintf(' At iteration %3.0f\n',k) X disp(' Partial solution (in ascending powers of s): ') X ddg = [1;ak]'; X ddg = fliplr(ddg); X for i = 1:r+1 X fprintf(' %18.14g',ddg(i)) X end X ksw = ksw + kdisp; X end X aprev = ak; X k = k + 1; X ak1 = [1;ak]; X for i = r:-2:1 X ak1(i) = - ak1(i); X end X v2 = - conv(ak1,ak1); X if ni > r+1 X v2 = [v2;zeros(ni-r-1,1)]; X end X v1 = zeros(1,ni-r); X v1(1) = v2(1); X mw = toeplitz(v2,v1); X m13 = mw(1:ni-r,:); X m23 = mw(ni-r+1:ni,:); X m33 = mw(ni+1:ni+r,:); X n13 = m13+h12*m23+h13*m33; X n23 = m23+h23*m33; X q = n13\e1; X ak = n21*(e2-n23*q); X% X% Determination of a candidate for the global minimum X% X bpom = -n31*ak - h33*m33*q; X jpo = energy(conv(num,[1,ak'])-conv(den,bpom'),conv(den,[1,ak'])); X if jpo > 0 X if jpo < jp X jp = jpo; X np = fliplr(bpom'); X dp = fliplr(ak'); X dp = [dp 1]; X kp = k; X end X end X delta = max(abs(ak-aprev)); Xend X% X% After a solution has been reached, the resulting polynomials coefficients X% are displayed according to ascending powers of s X% Xif ok == 1 X bk = -n31*ak -h33*m33*q; X nr = fliplr(bk'); X dr = fliplr(ak'); X dr = [dr 1]; X jr = jpo; X kr = k; X disp(' ') X fprintf(' Solution found in %3.0f iterations',k) X disp(' ') X disp(' Strike any key to continue ... '), pause X sour(nr,dr,jr,kr); X if kr ~= kp X disp(' A solution characterized by an index value lower than') X disp(' the one corresponding to the satisfaction of the ') X disp(' stopping criterion has been found') X disp(' Strike any key to continue ... '), pause X sour(np,dp,jp,kp); X end Xelse X disp(' ') X fprintf(' Cannot find any solution within %3.0f iterations',kmax) X disp(' ') X disp(' The best solution with respect to the index value follows') X disp(' Strike any key to continue ... '), pause X sour(np,dp,jp,kp); Xend END_OF_FILE if test 6809 -ne `wc -c <'l2siso.m'`; then echo shar: \"'l2siso.m'\" unpacked with wrong size! fi # end of 'l2siso.m' fi if test -f 'mino.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'mino.m'\" else echo shar: Extracting \"'mino.m'\" \(2608 characters\) sed "s/^X//" >'mino.m' <<'END_OF_FILE' Xfunction [num,den,mi,mo] = mino X% X% Description: X% Performs the input operations to get the values of the coefficients of X% the original model denominator and of the mo*mi matrix of the numerators X% (MIMO case). X% The order ni of the original model is set equal to the number of X% the denominator coefficients minus one. X% If the number y of the coefficients of some numerator polynomials X% specified by the user is less than ni, then this function adds further X% (ni-y) zero coefficients. X% X% Usage: X% [num,den,mi,mo] = mino; X% X% Output parameters: X% num - vector of the numerator coefficients stored by rows and in X% ascending order as follows X% [n(1,1,0) n(1,1,1) ... n(1,1,ni-1) ... n(1,2,0) ... n(1,2,ni-1) X% ... n(mo,mi,0) ... n(mo,mi,ni-1)] X% den - vector of the original denominator coefficients in ascending X% powers of s. If the denominator specified is not monic, then the X% function divides all the coefficients of the denominator and X% numerator by the coefficient d(ni) of the highest X% power of s. X% mi - number of inputs X% mo - number of outputs X% X Xsw = 0; Xwhile sw == 0 X disp(' ') X disp(' Specify the coefficients of the original model denominator ') X disp(' according to ascending powers of s, i.e. [d(0) d(1) ... ]. ') X den = input('? '); X% X% ni - order of the original model X% X ni = length(den) - 1; X if den(ni+1) == 0 X disp(' The coefficient of the term of highest degree is zero. Try again!') X else X sw = 1; X end Xend Xdisp(' ') Xfprintf(' The order of the original model is %3.0f\n',ni) Xmo = input('Specify the number of outputs. '); Xmi = input('Specify the number of inputs. '); Xdisp(' ') Xdisp(' Specify the numerator matrix row by row, giving the coefficients') Xdisp(' according to ascending powers of s, i.e. [n(i,j,0) n(i,j,1) ... ].') Xnum = []; Xfor i = 1:mo X disp(' ') X fprintf(' Row %3.0f\n',i) X for j = 1:mi X fprintf(' Column %3.0f',j) X sw = 0; X while sw == 0 X c = input('? '); X lc = length(c); X if lc > ni X fprintf ... X (' The coefficients must be less than or equal to %3.0f. Try again!',ni) X else X if lc < ni X cc = c; X c = zeros(1,ni); X c(1:lc) = cc; X end X num = [num,c]; X sw = 1; X end X end X end Xend X% X% Case in which den(s) is not monic X% Xif den(ni+1) ~= 1 X m = length(num); X x = den(ni+1); X for i = 1:ni+1 X den(i) = den(i)/x; X end X for i = 1:m X num(i) = num(i)/x; X end Xend END_OF_FILE if test 2608 -ne `wc -c <'mino.m'`; then echo shar: \"'mino.m'\" unpacked with wrong size! fi # end of 'mino.m' fi if test -f 'minr.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'minr.m'\" else echo shar: Extracting \"'minr.m'\" \(2586 characters\) sed "s/^X//" >'minr.m' <<'END_OF_FILE' Xfunction [ng0,dg0] = minr(mi,mo,ni) X% X% Description: X% Performs the input operations to get the order r of the reduced X% model and the first guesses for the denominator and numerators (MIMO case). X% If the number y of the coefficients of some polynomial of the numerator X% is less than r, then this function adds further (r-y) zero coefficients. X% X% Usage: X% [ng0,dg0] = minr(mi,mo,ni); X% X% Input parameters: X% mi - number of inputs X% mo - number of outputs X% ni - order of the original model X% X% Output parameters: X% ng0 - vector containing the first guess for the numerators stored by rows X% and in ascending order as follows X% [ng0(1,1,0) ng0(1,1,1) ... ng0(1,1,r-1) ... ng0(1,2,0) ... X% ng0(1,2,r-1) ... ng0(mo,mi,0) ... ng0(mo,mi,r-1)] X% dg0 - first guess for the reduced denominator. If this is not monic, X% then the function divides all its coefficients and all the X% coefficients of ng0 by dg0(r). X% X Xdisp(' ') X% X% r - order of the reduced model X% Xr = input('Specify the order of the reduced model. '); Xwhile r >= ni X fprintf(' The order must be less than %3.0f. Try again!',ni) X r = input('Specify the order of the reduced model. '); Xend X% X% First guess for the denominator X% Xsw = 0; Xwhile sw == 0 X disp(' Try a first guess for the denominator coefficients arranged in') X disp(' ascending order, i.e. [dg0(0) dg0(1) ... ]. ') X dg0 = input('? '); X if length(dg0) ~= r+1 X fprintf(' The length of this vector must be %3.0f. Try again!',r+1) X disp(' ') X elseif dg0(r+1) == 0 X fprintf(' The value of the term %3.0f is zero. Try again!',r+1) X disp(' ') X else X sw = 1; X end Xend X% X% First guess for the numerators X% Xdisp(' ') Xdisp(' Specify a first guess for numerators coefficients arranged in') Xdisp(' ascending order, i.e. [ng0(i,j,0) ng0(i,j,1) ... ].') Xng0 = []; Xfor i = 1:mo X disp(' ') X fprintf(' Row %3.0f\n',i) X for j = 1:mi X fprintf(' Column %3.0f',j) X sw = 0; X while sw == 0 X c = input('? '); X lc = length(c); X if lc > r X fprintf ... X (' The length must be less than or equal to %3.0f. Try again!',r) X else X if lc < r X cc = c; X c = zeros(1,r); X c(1:lc) = cc; X end X ng0 = [ng0,c]; X sw =1; X end X end X end Xend X% X% Case in which dg0 is not monic X% Xif dg0(r+1) ~= 1 X m = length(ng0); X x = dg0(r+1); X for i = 1:r+1 X dg0(i) = dg0(i)/x; X end X for i = 1:m X ng0(i) = ng0(i)/x; X end Xend END_OF_FILE if test 2586 -ne `wc -c <'minr.m'`; then echo shar: \"'minr.m'\" unpacked with wrong size! fi # end of 'minr.m' fi if test -f 'mour.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'mour.m'\" else echo shar: Extracting \"'mour.m'\" \(1384 characters\) sed "s/^X//" >'mour.m' <<'END_OF_FILE' Xfunction mour(n,d,j,k,mi,mo) X% X% Description: X% Performs the output operations for displaying the denominator, X% the numerator, the value of the index and the iteration number X% (MIMO case). X% X% Usage: X% mour(n,d,j,k,mi,mo); X% X% Input parameters: X% n - reduced order numerator X% d - reduced order denominator X% j - index value X% k - iteration number X% mi - number of inputs X% mo - number of outputs X% X Xdisp(' ') X% X% r - order of the reduced model X% Xr = length(d) - 1; Xfprintf(' The order of the reduced model is %3.0g\n',r) X% X% Iteration number X% Xfprintf(' Iteration number %3.0f\n',k) X% X% Index value X% Xfprintf(' The index value is equal to %20.14g\n\n',j) X% X% Denominator coefficients X% Xdisp(' Denominator in ascending powers of s. ') Xfor i = 1:r+1 X fprintf(' %18.14g',d(i)) Xend X% Xdisp(' ') Xdisp(' Strike any key to continue ... '), pause X% X% Numerator coefficients X% Xdisp(' ') Xdisp(' Numerator in ascending powers of s. ') Xfor i = 1:mo X disp(' ') X fprintf(' Row %3.0f\n',i) X st = (i-1)*mi*r; X for j = 1:mi X fprintf(' Column %3.0f\n',j) X li = st + (j-1)*r + 1; X ui = li + r - 1; X for l = li:ui X fprintf(' %18.14g',n(l)) X end X disp(' ') X if i == mo & j == mi X disp(' Strike any key to stop ... '), pause X else X disp(' Strike any key to continue ... '), pause X end X end Xend END_OF_FILE if test 1384 -ne `wc -c <'mour.m'`; then echo shar: \"'mour.m'\" unpacked with wrong size! fi # end of 'mour.m' fi if test -f 'rever.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'rever.m'\" else echo shar: Extracting \"'rever.m'\" \(1827 characters\) sed "s/^X//" >'rever.m' <<'END_OF_FILE' Xfunction [nrev,drev] = rever(num,den,mi,mo,i,j) X% X% Description: X% It reverses the order of the denominator and of the numerators X% coefficients (from ascending to descending powers of s). X% If the user supplies mi,mo,i and j (input parameters), then the X% output parameter nrev is a vector that contains all the coefficients X% for the polynomial corresponding to the (i,j)-th position in the X% numerator matrix (in descending powers of s). X% In the MIMO case, if nrev only is specified as output parameter, X% then the numerator only is reversed. X% X% Usage: X% [nrev,drev] = rever(num,den,mi,mo,i,j); X% X% Input parameters: X% num - numerators coefficients in ascending powers of s X% den - denominator coefficients in ascending powers of s X% Optionals: X% mi - number of inputs X% mo - number of outputs X% i,j - position of the polynomial entry in the numerator matrix X% X% Output parameters: X% nrev - numerator coefficients in descending powers of s, stored X% according to the above Description X% drev - denominator coefficients in descending powers of s X% X Xdisp(' ') X% X% n - order of the model X% Xn = length(den) - 1; Xnl = length(num); X% X% Check whether SISO or MIMO X% Xif nargin < 3 X% X% The coefficients order of the denominator and of the numerator X% is reversed (SISO case) X% X drev = fliplr(den); X nrev = fliplr(num); Xelseif nl ~= mi*mo*n X fprintf (' The length of the numerator is not %3.0f\n',mi*mo*n) X disp(' Controls the values mi and mo') Xelseif i > mo & j > mi X disp(' Controls the values i and j') Xelse X if nargout == 2 X% X% The coefficients order of the denominator is reversed (MIMO case) X% X drev = fliplr(den); X end X% X% The (i,j)-th numerator is stored (MIMO case) X% X ist = (i-1)*mi*n + (j-1)*n; X for i = 1:n X nrev(n-i+1) = num(ist+i); X end Xend END_OF_FILE if test 1827 -ne `wc -c <'rever.m'`; then echo shar: \"'rever.m'\" unpacked with wrong size! fi # end of 'rever.m' fi if test -f 'sino.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'sino.m'\" else echo shar: Extracting \"'sino.m'\" \(1908 characters\) sed "s/^X//" >'sino.m' <<'END_OF_FILE' Xfunction [num,den] = sino X% X% Description: X% Performs the input operations to get the values of the original X% model denominator and numerator (SISO case). X% The order ni of the original model is set equal to the number of X% the denominator coefficients minus one. X% If the number y of the numerator coefficients specified by the user X% is less than ni, then this function adds further (ni-y) zero X% coefficients. X% X% Usage: X% [num,den] = sino; X% X% Output parameters: X% num - original model numerator (ascending powers of s), X% den - original model denominator (ascending powers of s). If the X% denominator specified is not monic, then the function divides X% all the coefficients of the denominator and numerator by the X% coefficient d(ni) of the highest power of s. X% X Xsw = 0; Xwhile sw == 0 X disp(' ') X disp(' Specify the coefficients of the original model denominator') X disp(' according to ascending powers of s, i.e. [d(0) d(1) ... ]. ') X den = input('? '); X% X% ni - order of the original model X% X ni = length(den) - 1; X if den(ni+1) == 0 X disp(' The coefficient of the term of highest degree is zero. Try again!') X else X sw = 1; X end Xend Xdisp(' ') Xfprintf(' The order of the original model is %3.0f\n',ni) Xdisp(' ') Xdisp(' Specify the coefficients of the original model numerator') Xdisp(' according to ascending powers of s, i.e. [n(0) n(1) ... ]. ') Xsw = 0; Xwhile sw == 0 X num = input('? '); X lc = length(num); X if lc > ni X fprintf ... X (' The coefficients must be less than or equal to %3.0f. Try again!',ni) X else X if lc < ni X cc = num; X num = zeros(1,ni); X num(1:lc) = cc; X end X sw = 1; X end Xend X% X% Case in which den(s) is not monic X% Xif den(ni+1) ~= 1 X x = den(ni+1); X for i = 1:ni+1 X den(i) = den(i)/x; X end X for i = 1:ni X num(i) = num(i)/x; X end Xend END_OF_FILE if test 1908 -ne `wc -c <'sino.m'`; then echo shar: \"'sino.m'\" unpacked with wrong size! fi # end of 'sino.m' fi if test -f 'sinr.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'sinr.m'\" else echo shar: Extracting \"'sinr.m'\" \(1266 characters\) sed "s/^X//" >'sinr.m' <<'END_OF_FILE' Xfunction dg0 = sinr(ni) X% X% Description: X% Performs the input operations to get the order r of the reduced X% model and the first guess for its denominator (SISO case). X% X% Usage: X% dg0 = sinr(ni); X% X% Input parameter: X% ni - order of the original model X% X% Output parameter: X% dg0 - first guess for the reduced denominator. If this is not monic, X% then the function divides all its coefficients by dg0(r). X% X Xdisp(' ') X% X% r - order of the reduced model X% Xr = input('Specify the order of the reduced model. '); Xwhile r >= ni X fprintf(' The order must be less than %3.0f. Try again!',ni) X r = input('Specify the order of the reduced model. '); Xend X% X% First guess for the denominator X% Xsw = 0; Xwhile sw == 0 X disp(' Try a first guess for the denominator coefficients arranged in ') X disp(' ascending order, i.e. [dg0(0) dg0(1) ... ]. ') X dg0 = input('? '); X if length(dg0) ~= r+1 X fprintf(' The length of this vector must be %3.0f. Try again!',r+1) X disp(' ') X elseif dg0(r+1) == 0 X fprintf(' The value of the term %3.0f is zero. Try again!',r+1) X disp(' ') X else X sw = 1; X end Xend X% X% Case in which dg0 is not monic X% Xif dg0(r+1) ~= 1 X x = dg0(r+1); X for i = 1:r+1 X dg0(i) = dg0(i)/x; X end Xend END_OF_FILE if test 1266 -ne `wc -c <'sinr.m'`; then echo shar: \"'sinr.m'\" unpacked with wrong size! fi # end of 'sinr.m' fi if test -f 'sour.m' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'sour.m'\" else echo shar: Extracting \"'sour.m'\" \(988 characters\) sed "s/^X//" >'sour.m' <<'END_OF_FILE' Xfunction sour(n,d,j,k) X% X% Description: X% Performs the output operations for displaying the denominator, X% the numerator, the value of the index and the iteration number X% (SISO case). X% X% Usage: X% sour(n,d,j,k); X% X% Input parameters: X% n - reduced order numerator X% d - reduced order denominator X% j - index value X% k - iteration number X% X Xdisp(' ') X% X% r - order of the reduced model X% Xr = length(d) - 1; Xfprintf(' The order of the reduced model is %3.0g\n',r) X% X% Iteration number X% Xfprintf(' Iteration number %3.0f\n',k) X% X% Index value X% Xfprintf(' The index value is equal to %20.14g\n\n',j) X% X% Denominator coefficients X% Xdisp(' Denominator in ascending powers of s. ') Xfor i = 1:r+1 X fprintf(' %18.14g',d(i)) Xend X% Xdisp(' ') Xdisp(' Strike any key to continue ... '), pause X% X% Numerator coefficients X% Xdisp(' ') Xdisp(' Numerator in ascending powers of s. ') Xfor i = 1:r X fprintf(' %18.14g',n(i)) Xend Xdisp(' ') Xdisp(' Strike any key to stop ... '), pause END_OF_FILE if test 988 -ne `wc -c <'sour.m'`; then echo shar: \"'sour.m'\" unpacked with wrong size! fi # end of 'sour.m' fi echo shar: End of shell archive. exit 0