MARSYAS is a computer-aided control system design package for the simulation and analysis of dynamic systems. In the summer of 1991 MARSYAS was updated to allow for the analysis of sampled-data systems in terms of frequency response, stability, etc. This update was continued during the summer of 1992 in order to extend further MARSYAS commands to the study of sampled-data systems. Further work was done to examine the computation of OPENAT transfer functions, root-locii and w-plane frequency response plots

Hodel, A. Scottedward

English

NASA Technical Reports Server

N9 ,,
1992
NASA/ASEE SUMMER FACULTY FELLOWSHIP PROGRAM
Numerical Methods for the Analysis of Sampled-Data Systems and for the Computation of
System Zeros
Prepared by:
Academic Rank:
Institution
NASA/MSFC:
Office:
Division:
Branch:
MSFC Colleague:
Contract No:
A. Scottedward Hodel, Ph. D.
Assistant Professor
Department of Electrical Engineering
Auburn University
Structures and Dynamics Laboratory
Control System Division
Mechanical Systems Control Branch
D. P. Vallely
NGT-01-008-021
The University of Alabama at Huntsville
XXI
https://ntrs.nasa.gov/search.jsp?R=19930008111 2020-03-17T09:11:07+00:00Z

1 Introduction
MARSYAS is a computer-aided control system design package for the simulation and analysis of dynamic
systems. In the summer of 1991 MARSYAS was updated to allow for the analysis of sampled-data systems
in terms of frequency response, stability, etc. This update was continued during the summer of 1992 in order
to extend further MARSYAS commands to the study of sampled-data systems. Further work was done to
examine the computation of 0PENAT transfer functions, root-locii and w-plane frequency response plots.
2 Sampled-data systems with feed-forward coefficients
Consider a sampled data system whose state-space equations are
z¢(t) = Acx_(t) + Bodyd(kT) + B_u(t)
y(t) = Ccxc(t) + Dcuyd(kT) + Dcu(t)
,,u(t) = cuo o(t) + Dduyu(kT)+ Duou(t)
(2.1)
(2.2)
(2.3)
with corresponding discrete time subsystem
zd(kT + T) = Adzd(kT)+ Buud(kT) (2.4)
yd(kT) = Cdzd(kT) + Ddud(kT). (2.5)
(Purely discrete time inputs ra(kT) to the discrete time system may be incorporated into the above equations
by augmenting the input vector u(t) and the feed-forward matrix Ddc.) The feed-forward coefficient Dd
causes the states z¢(t) to depend not only on their continuous values through the coefficient matrix A¢, but
also on their sampled values through the coefficient matrix BcaDdCd_. As of MARSYAS version 6.0.5, this
dependence is dealt with by creating a matrix ZX = BcdDdCdc in the linearized continuous time dynamics.
Unfortunately, MARSYAS version 6.0.5 does not correctly treat the case where there are feed-forward
coefficients in both the continuous time and discrete time blocks (i.e. Dad ¢ 0 and Dd # 0). In this case,
MARSYAS infers an algebraic loop where there is none. For example, consider the sampled-data system
with continuous time block
y(t) = yu(kT) ud(t) = u(t)/2- yd(kT)
and discrete time block
yd(kT) = kdud(kT).
The code has no (explicit) states, yet because of the direct feed-forward gains (i.e., the "D" matrices) in the
discrete and continuous time blocks of the linear system, a "hidden state" associated with the A/D converters
in the system manifests itself at sampling times. The MARSYAS-6.0.5-generated system equations are
W [ 1] = (2.S00000Z-01) * U [ 1]
while the simulated output values display "spikes" at sampling times and the values between sampling times
incorrect. There is a specific order in which the simulation updates must occur at each "scheduled" run of the
discrete time systems: (1) Evaluate the continuous time system states and algebraic outputs. (2) Update the
discrete time system states and algebraic outputs. (3) the algebraic outputs of the continuous time system to
reflect the discrete time system changes. Once these three steps are completed, the simulation may proceed
with the next integration step. Similarly, the analysis of discrete time systems must be modified as follows.
Lemma 2.6 Let a sampled data system be defined by the equations (_.I)-(_.5) with sampling time T given.
Then
ud(kT) = Cdczc(kT) + Detcu(kT) + Detdyd(kT - T);
that is, the discrete-time output yet becomes a system state when Ddd # O.
XXI-1
Figure1: OPENATlinearized system
Based on the above operation, the vector y_(t) can be used as a vector state that replaces the "ZX" input
used in MARSYAS version 6.0.5 as follows. Let z_y(kT) = yd(kT -- T). Then the overall system simplifies
to
xc(t) = Acz¢(t) + Bcdyd(kT) + B¢u(t)
y(t) = Cezc(t) + Dcayd(kT) + Dcu(t)
u_(t) = cdoxo(t) + D_¢u(t)
x_(kT + T) : Adr_(kT) + B_ud(kT) + BdDdaxdy(kT)
xdy(kT + T) = Cdczd(kT) + D_Dddxdy(kT) + Daud(kT)
yd(kT) = C_z_(kT) + D_Dd_x_u + D_ud(kT).
3 OPENAT analysis: options and an example
The 0P_.liAT command allows MARSYAS users to examine system robustness with respect to an individual
parameter by breaking a specified signal path in a closed-loop system model and examining the poles/zeros
of the newly created open loop system.
The current MARSYAS implementation (and the original MARSYAS implementation) treat the 0P_-NAT
command as a special case of the IN0trr command. The manual for the original MARSYAS implementation
describes an alternate technique for computing the 0PENAT transfer function; this alternate technique forms
the basis of the root locus calculation, and is thus of interest in both of these calculations. The two techniques
are described below. A simple numerical example is given that indicates that algebraic loops (i.e., a non-zero
D matrix in a linearized OPENAT system) can render the alternate algorithm inaccurate; it is advised that
the Ill0trr approach be used in both 0P_.NAT and root locus calculations.
The 0P_.liAT command can be summarized as follows. Given a closed loop continuous time system
= f(x, w) with algebraic constraint 0 = g(x, w), 0PEIqAT k selects a gain value in the system and linearizes
about an operating point in order to obtain the SISO system shown in Figure 1. The 0P_.IlAT command
identifies the poles and zeros of the transfer function n(s)/d(s). The IN0trr algorithm computes the poles
and zeros by computing the matrices (A, B, C, D) that characterized a state-space realization of n(s)/d(s)
and then computing the (finite)generalized eigenvalues of the matrix pencil ([-A-B ] [I 0])C D -A 0 "
The alternate closed-loop matrix method is as follows For a fixed linearization, let A(k) be the system
Jacobian matrix evaluated with K = k; i.e. A(k) = (A + k/(1 - Dk)BC). Then det(sI - A(k)) =
a(k)(d(s) - kn(s)) where a(k)-1 is the leading coefficient of (d(s) - kn(s)). Clearly, the poles of the 0PEllAT
system are the eigenvalues of A(0), since det(sI - A(0)) = a(O)(d(s) - O. n(s)) = d(s).
In order to obtain the zeros of the system, it is required to compute the algebraic gain D and the closed
loop matrix A(1) = (A(0) + 1/(1 - D)BC). Let/}, C be a column and row vector, respectively, such that
/}C = (D - 1)(A(1) - A(0); then the zeros of the 0PENAT transfer function may be obtained by computing
XXI-2
INOUTmethod Closed-loop matrix method
I I I Ib -Dli
1 _l.Oe-13 1.9984e-14 1 l.Oe-13 2"1.9984e-14
0.I000 ,_,l.Oe-14 1.1102e-14 0.I000 _l.Oe-12 2"I.II02e-14
0.0100 _1.0e-12 9.9920e-15 0.0100 _l.0e-lO 2"9.9920e-15
1.0000e-03 _,1.0e-12 9.9920e-15 l.O000e-03 _1.0e-09 2"9.9920e-15
1.0000e-04 ,_I.0e-ll 9.9920e-15 1.0000e-04 _,1.0e-07 2"9.9920e-15
1.0000e-05 _1.0e-ll 9.9920e-15 1.0000e-05 _1.0e-06 2"9.9920e-15
1.0000e-06 _1.0e-10 9.9920e-15 1.0000e-06 _1.0e-04 2"9.9920e-15
1.0000e-07 _1.0e-09 9.9920e-15 1.0000e-07 _,1.0e-03 2"9.9920e-15
Figure 2: Resultsoftestsofopenat algorithms
the finite generalized eigenvalues of the matrix pencil C D 0 0 "
matrix method discussed in the manual for the original MARSYAS implementation.
The discussion of the closed-loop matrix method in the original MARSYAS manual did not correctly
treat the case of D ¢ 0. If D ¢ 0, in particular, if D _ 1, then the Ii0trr method is clearly superior to the
closed-loop matrix method, as shown in the following example.
Example 3.1 Let the OPENAT linearized system have coefficient matrices A =
C= [ 2 1 ],andD=l+e. The exact system zeros can be found to be
± 1+y4- -
-I -I ' 2 '
Observe that as e --_ 0, the zeros coincide at A = 3.
A computed linearization (.4, ]_, C, D) = (A, B, C, D) will typically have roundoff noise in each matrix in
the matrix norm; i.e., ,,[i_ -All < tt IIAI] where tt is a function of machine precision and theproportion to
algorithm condition. The example was selected so that matrix entries were "balanced;" that is, the matrix
entries were of approximately the same magnitude so that round-off effects on problem condition would be
reduced. For the purposes of a numerical experiment, small error (_ 10 -14) was deliberately introduced into
the problem in order to examine algorithm sensitivity. Since the computed value of Ab is near to the correct
computed value, this method yields good results.
The closed-loop matrix method involves computing A(0), A(1), and the open-loop D-matrix. it was not
assumed that the gain D involved in computing A(1) = A-1/(D - 1)BC is exactly the same as that gained
from an IIOOT linearization; i.e., IHolrr computes (A, B, C, D1) - (A, B, C, D), while A(1) was computed as
A(1) = A + 1/(/9_ - 1)BC. The respective results are shown in Figure 2. It is clear that the IiOtrr method
yields superior accuracy; it is important to recognize that the difference in these two methods is a result of
the difference between 1/(D1 - 1) and 1/(D) - 1); since the estimated values of/} and C in the closed loop
matrix method will be scaled by the ratio (DI - 1)/(192 - 1) from their "correct" values, it is to be expected
that the closed-loop matrix method will yield poor results when D _ 1.
E]
On the basis of the above example and analysis, it is recommended that both 0PENAT and root locus
calculations be performed on the basis of the matrices (A, B, C, D) obtained from an IlIotrr linearization of
the appropriate system. FORTRAN code for the mvzero routine and the associated numerical balancing
procedure [2] have heed delivered to John Tiller, BCSS, this summer for use in MARSYAS.
XXI-3
The calculationof OPE]IATtransferfunctionsand root locus data become much more complicated in a
sampled-data system. Ifthe OPENAT operation isperformed on a gain k in the discretetime portion of the
system, then the OPEIIATfunction may be executed in the normal fashionon the equivalentdiscretetime
system. However, ifthe OPENAT isperformed on a gain k inthe continuoustime part ofthe system, then the
OPENAT function ceasesto have a clearmeaning. Ifan artificialcontinuous time input u and output y are
placed around the gain,then there isno rationaltransferfunctionfrom u to y because of the time-varying
dependence on the discretetime subsystem. Instead,discretetime input-outputpair isinsertedabout the
continuous time gain element K.
This system is not amenable to an OPEIAT style analysis, since the transfer function varies in a transcen-
dental (non-algebraic) fashion with the gain value K; closed form perturbation analysis of the equivalent
discrete time system is intractable since the matrices A and BC do not commute in general. (eAsc =
eAe Bc = eBCe A _ A(BC) = (BC)A.) Hence a single "open-loop" 0PENAT transfer function calculated
with K = 0 is not meaningful to the user. However, a "transcendental root-locus" from u to y may be
computed by computing the poles and zeros of an equivalent discrete-time plant for various values of k.
Closed-form analysis of transcendental root-iocii is hindered by the property that e(Ac+KBcD°)t cannot
be computed from eA_t and e(Bcc¢)t except when Ac commutes with BcCc; this is clearly not the case
in practice. It may be possible to gain some norm-bounds on the closed-loop eigenvalues of the overall
discrete-time system in terms of the feedback gain K, but this is an open question.
4 W-plane analysis
MARSYAS currentlyprovides z-planefrequency response plotsof discrete-timesystems, w-plane analysis
isoften used in practice to allow designersto employ continuous time design techniques to discrete-time
systems; see, e.g., [1]. G(z) --_ G(w) by w = _ _ . s-plane, z-plane, and w-plane frequency response
2 e'r-I 2 tanh (_) So as s = jw follows the imaginary axis,plots are related by z -- e aT and w = Te,-rrW-T+l= Iv
z = ej'°T follows the unit circle (periodic with period 2_r/T), and
w = _ tanh = -_ tan
follows the entire imaginary axis each time z = ej"T rotates about the unit circle The inverse bilinear
1 + wT/2
transform from w to z is z = 1 -wT/2; hence the Nyquist plots in the z and w plane will be identical.
Magnitude and phase information in the z and w-planes will he the same except for a "warping" of the
zx tan(wT/2).w-plane frequency variable u =
References
[1] G. F. Franklin, J. D. Powell, and M. L. Workman. Digital Control of Dynamic Systems. Addison-Wesley,
2nd edition, 1990.
[2] A. Scottedward Hodel and John Tiller. Computation of system zeros with balancing. In Proceedings of
the 29 th Allerton Conference on Comm., Contr., and Computing, Monticello, II., Oct 2-4 1991.
XXI-4


Numerical methods for the analysis of sampled-data systems and for the computation of system zeros

Abstract

Similar works

Full text

Available Versions

NASA Technical Reports Server