#! /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 <file", e.g.. If this archive is complete, you
# will see the following message at the end:
# "End of shell archive."
# Contents: READ.ME energy.m l2mimo.m l2siso.m mino.m minr.m mour.m
# rever.m sino.m sinr.m sour.m
# Wrapped by michela@coco on Wed Sep 27 12:14:48 1995
PATH=/bin:/usr/bin:/usr/ucb ; export PATH
if test -f 'READ.ME' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'READ.ME'\"
else
echo shar: Extracting \"'READ.ME'\" \(1820 characters\)
sed "s/^X//" >'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