.4 ”$.33

«.4
am.“

.53 ~

 

1):...

l
‘5
h

.3
h
a
3
1

12!

M“.

I . .«Wu‘ifimw
.4. flu. u.
'OI‘A I
i}: .. A!
rdusaun.

{14

a“
I .stuxfiuflv
, 3F “Jaguxv
JFK); \fitflfin‘

I .V f0
5 r! . 5.5.4.
a... (3va . 0.0.».
.3: . .

23
51v

 

{Mr

M

.I .. .v.. . ..
2.. . flan. 33‘

mm. fie? 1 . 3§i§§§

 

This is to certify that the
dissertation entitled

NETWORKED ASSEMBLY OF NONLINEAR PHYSICAL
SYSTEM MODELS

presented by

Elliot Motato

has been accepted towards fulfillment
of the requirements for the

PhD. degree in Mechanical Engineering

 

 

_ Major Professor Milyturesfi7
Dal te

MSU is an Affirmative Action/Equal Opportunity Institution

 

~i- .—.-.—-—-—.--«-vup.-.-.—--.-.---.-.-—-

 

 

UBRARY
Michigan State
University

 

 

 

PLACE IN RETURN BOX to remove this checkout from your record.
TO AVOID FINES return on or before date due.
MAY BE RECALLED with earlier due date if requested.

 

DATE DUE DATE DUE DATE DUE

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

6/07 p:/ClRC/DateDue.indd-p.1

 

 

NETWORKED ASSEMBLY OF NONLINEAR PHYSICAL
SYSTEM MODELS

By

Elliot Motato

A DISSERTATION

Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

DOCTOR OF PHILOSOPHY

Department of Mechanical Engineering

2007

ABSTRACT

NETWORKED ASSEMBLY OF NONLINEAR PHYSICAL
SYSTEM MODELS

By

Elliot Motato

Global engineering design requires efficient model integration. This process is char-
acterized by the capability of performing an automatic and recursive model assembly.
Model integration is a complicated issue when complex assemblies are involved. Model
reuse is the key for handling assembly complexity. The ability to reuse and exchange
models relies on a standard format [1]. This condition is fundamental to avoid model
reformulation.

The Modular Modeling Method (MMM) [2] is a recursive model integration al-
gorithm designed to achieve automatic integration of linear models. MMM linear
standard models are dynamic matrices with a unique input-output variable represen—
tation. MMM characteristics facilitates the networked integration of large complex
linear models.

The objective of this work is to extend the model assembly properties of MMM to
assemble nonlinear physical models. This work can divided in two parts. In the first
part, the problem of assembling multi-port nonlinear physical models that perform
at a constant operating point is solved. In the second part, the problem of assem-
bly nonlinear physical models that perform at a region of operation is solved. Once
achieved, the MMM nonlinear methodology will contribute to enable cooperative
global engineering design, by facilitating the automatic integration and the simula-
tions of nonlinear physical engineered artifacts which are designed and provided from

different and possibly geographical distant locations.

To my dears Sandra and Sarah.

iii

ACKNOWLEDGMENTS

I would like to thank my Advisor Dr. Clark Radcliffe for his continuous guidance

and endless support during my years at Michigan State University.
My thanks to Dr. Khalil, Dr. Gokcek, Dr. Sticklen and Dr. Mukerjee for serving
as members of my PhD committee. They provide useful suggestions to improve this

work.

My thanks to Umar F arooq, Jimmy Issa, Jeff Roads, Nanda Methil, Jorge Villa and

Wesley Zanardelli for providing me valuable help through their friendship.

Finally, I would like to thank the graduate secretary Aida Montalvo for her help.

iv

TABLE OF CONTENTS

Introduction ................................ l
1.1 Global Engineering Design ........................ 1
1.2 Contributions of This Work ....................... 4
1.3 Limitations of the Nonlinear l\i'Il\/Il\-/I Extension ............. 7
Background ................................ 9
2.1 Internet Distribution of Models ..................... 9
2.2 The Standard Linear Physical Model Formats ............. 11
2.3 Physical Assembly Constraints ...................... 14
2.4 The Modular Modeling Assembly Procedure .............. 20
2.5 Example of 3 Linear MMM Assembly .................. 24
2.6 Summary ............ ’ ..................... 2 8
Assembly of Affine Physical System Models ............. 30
3.1 Affine Physical System Models ...................... 30
3.2 Internet Distribution of Affine Models .................. 36
3.3 The Standard Model Formats ...................... 37
3.4 Physical Assembly Constraints ...................... 39
3.5 Assembly of Affine Physical Models ................... 40
3.6 Example .................................. 42
3.7 Summary ................................. 47
Volterra Model Formats ......................... 49
4.1 Background ................................ 49
4.2 Port Based Nonlinear ODE Models ................... 50
4.3 Volterra Models .............................. 51
4.4 The Nonlinear Multiplicative Operator ................. 57
4.5 Generating Volterra Models From Port-Based ODEs .......... 60
4.6 Example .................................. 67
4.7 Summary ................................. 73
Assembly of Nonlinear Physical Systems ............... 76
5.1 The MMM Algorithm to Assemble Nonlinear Physical Models . . . . 76
5.2 Simulation of the Assembled Nonlinear Physical Model ........ 78
5.3 Example .................................. 79
5.4 Summary ................................. 89
Conclusions ................................ 90

APPENDICES ................................ 92

A Singular Internal Stiffness Matrix ................... 92
B The Assembly Operator ......................... 94
BIBLIOGRAPHY .............................. 96

vi

CHAPTER 1

Introduction

1.1 Global Engineering Design

Engineering analysts are in the process of creating a global engineering strategy that
is able to integrate product design, product development, marketing analysis, and
manufacturing process [3]. Today’s engineers can communicate with of suppliers
through supply chain management systems over the Internet [4]. The trend is to use
the Internet as the media to achieve this global engineering design strategy.

Product design in the context of a global engineering strategy demands the assem-
bly of dynamic system models from component models retrieved over the Internet.
This process requires four important characteristics: 1. Multi-energy domain compo-
nent models must have a unique standard format, 2. The exchange of model informa-
tion must be executed in a single-query transmission, 3. The models must describe
only external behavior, and 4. The assembly process must be recursive. These four
characteristics make global engineering design practical.

A unique, standard model format is the key to handling model exchange through
model reuse [4]. A modeling methodology using a unique model representation facil-
itates model query standardization. A unique model format query prevents model
reformulation and decreases model exchange time computation. The Finite Element

Method (FEM) [5], is a modeling methodology that uses a unique model format. FEM

uses the same format to represent both elements and the system assembled from those
elements. A modeling method that uses a unique standard format is called modular.

Single—query exchange of model information reduces network traffic during the
assembly of dynamic models through the Internet. A single-query exchange retrieves
the full component model using a single request and answer on the Internet. Repeated
queries increase network traffic dramatically. The Differential Algebraic Equation
(DAE) approach [6] is an example of a. model format, which can be exchanged in a
single—query network transmission as a set of ordinary differential equations (ODEs)
and algebraic constraints. Global engineering strategy demands assembly processes
that reduce network load through exchange of model information with a single-query.
A model assembly method that has this characteristic is called single-query.

An input-output model predicts external behavior using only input and output
variables to protect internal proprietary design details. Because design is a dominant
cost of new product development, internal product design details must be protected
from competitors. These design details might include, the components used in the
assembly, the order of connection of those components, the physical parameters of l
the components and the performance of each component. Protection of pr0prietary
information is critical to the commercial acceptance of any model exchange system.
Gu and Asada in [4] have used input and output variables to allow the co—simulation
of a collection of dynamic sub-simulators without disclosing proprietary information.
A model assembly method that uses models that predict external behavior using
external input and output variables is called esrtemal.

A recursive model assembly process uses standard format component models to
produce an assembly model in the same format. Once model assembly recursion
is established at a single level, the model assembly method is easily extended to
higher-level, more complex system models. DAB [6] is an example of an assembly

methodology that recursively obtains DAE system models from either DAE elements

or DAE subsystems. A model assembly process that uses standard format component
models to produce an assembly model in the same format is called recursive.

All four characteristics are simultaneously required for a successful global engineer-
ing model assembly method. Co—simulation, [4], is external but is not single-query
because network iteration is required to run dynamic sub-simulators. The DAE ap-
proach in [6], is recursive and single-query but is not external because DAE models
provide internal information about the assembly components, component connectiv-
ity and internal parameters. Bond graphs [7] are modular and recursive but not
external because they provide information about assembly components, component
connectivity and internal parameters. FEM [5] is modular and single-query but not
recursive. Assembly of FEM models generally requires global reformulation to guar-
antee geometric nodal compatibility. A multi-energy domain, modular, single-query,
external and recursive, model assembly methodology is needed.

MMM [8]-[9] is modeling strategy that satisfies global engineering design require-
ments. MMM uses two standard model formats. Dynamic matrices are used in the
networked distribution and assembly of linear models. Transfer function matrices are
used in model simulation. MMM is characterized by the advantages that result from
using these two model representations. This work is an extension of the MMM for
assembly of nonlinear physical models represented through Volterra models. Using
this extension, nonlinear dynamic models of assemblies can be built and distributed
while hiding the topology and characteristics of their structural subassemblies.

This work is divided in four parts. In the first part, the linear MMM procedure,
the standard linear model formats and the physical constraints to assemble physi-
cal models are reviewed. In the second part, the problem of assembling nonlinear
physical models that operate at a constant point is solved. In the third part, a pro-
cedure to obtain Volterra models from port based nonlinear differential equation is

explained. Finally in the fourth part, the MMM procedure for assembling Volterra

models performing at a region of operation is derived.

1.2 Contributions of This Work

This work makes six (6) significant, original contributions to mechanical engineering.

These contributions are:

1. An extension of the Modular Modeling Methodology (MMM) to nonlinear systems.
2. The recognition that two model formats are necessary for physical system model
assembly and simulation.

3. A method to assemble physical port-based affine ODEs around an equilibrium
operating point.

4. A method to obtain subsystems operating points from the system assembly oper-
ating point outputs.

5. A method to obtain frequency domain Volterra system models from MIMO port-
based ODEs.

6. A method to assemble physical MIMO port—based Volterra system models.

Each of these contributions is individually important to the development of a Global
Engineering Design (GED) strategy.

The extension of the MMM to nonlinear systems is important. This extension
allows the application of the MMM to assemble, distribute and simulate system mod-
els for a large class of nonlinear physical systems. This class is consists of nonlinear
physical systems that are analytic or that have fading memory. Before this extension
was developed, the application of the MMM was limited to the assembly, distribu-
tion and simulation of linear physical system models. The extension of the MMM to
nonlinear systems is discussed throughout the thesis. Model distribution is discussed
in the Chapter Two. Nonlinear subsystem models assembly is discussed in Chapters

Three, Four and Five. Finally, nonlinear assembled system simulation is discussed at

the end of the Chapter Five.

The recognition that two model formats are necessary for physical system model
assembly and simulation is important. As discussed below, the use of two model
formats decreases the computational time required to assemble and simulate physical
system models. These two model formats are the linear/ nonlinear transfer function
representation and the linear/ nonlinear dynamic representation. The transfer func—
tion representation is the traditional format for simulation. This thesis shows for the
first time, that in a constraint-based model assembly, the Transfer Function represen-
tation is not appropriate for MMM assembly because a matrix inversion is required.
Matrix inversion demands substantial computational time, a clear disadvantage in a
global distributed environment. Additionally, for this inverse to exist, system models
are limited to have independent inputs and independent outputs. This is a big limi-
tation on allowable models that excludes many dynamic system models. In contrast,
if the format used for assembly is the dynamic representation, no matrix inversion is
required and system models with dependent inputs and dependent outputs can also
be assembled. The recognition that the two model formats are necessary for efficient
physical system model assembly and simulation is discussed in the Chapter 2.

A method to assemble physical port-based affine ODEs around an equilibrium op-
erating point is important. Affine systems often result of local linearization about an
operating point. In this case, the local system model is linear in deviation variables
and non-linear in physical variables. Because assembly constraints are always given
in terms of physical variables, an explicit method to apply physical assembly con-
straints for affine physical models was first developed in this work. This new method
addresses one of the most common nonlinear systems in mechanical engineering. Us-
ing this method many practical physical nonlinear system models can be obtained
by assembling the affine approximations of their subsystems models. In addition,

this method is computationally efficient because it is recursive and does not require

matrix inversions. These characteristics make this affine assembly procedure original
and ideal for being used in the MMM. A method to assemble physical port-based
affine ODEs around an equilibrium operating point is discussed in the Chapter 3.

The explicit method to obtain subsystems operating points from the system assem-
bly operating point outputs is important. This method for the first time provides an
explicit, closed form solution to the general operating point problem. In a recursive,
closed, form, the method specifies information required to roll-down an operating
point specification through every subsystem to all the lowest level components of a
system. For each lowest level component, the method provides the exact port output
values where the operating point representation of the nonlinear component model
should be developed. The method then provides an explicit method for computing
operating point inputs of the component level and use these inputs to compute the
operating point inputs of the assembled system model. This method is demonstrated
through Affine approximations, one of the standard formats used in the MMM nonlin-
ear extensions to nonlinear system models. Determining the operating point inputs
and outputs of nonlinear subsystems in a system assembly is a classical problem [4].
This method provides an explicit solution to this problem for an important class of
mechanical engineering system models: port-based physical system assemblies. The
method to obtain subsystems Operating points from the operating outputs of an as-
sembly is discussed in Chapter 3.

The method developed here to generate frequency domain Volterra system models
from MIMO external port-based ODEs is important. Even though a similar procedure
for the SISO case is available in the literature [10], a method for MIMO port-based,
external nonlinear systems has not been published. The MIMO procedure requires
a nonlinear vector operator to obtain the frequency domain kernels. This nonlinear
operator is not required in the SISO case. An advantage of using this MIMO method

is that the two standard MMM nonlinear formats can be easily obtained. The first

format is the Volterra transfer function representation. This format is traditionally
used in model analysis and simulation. The second format is the Volterra dynamic
representation. This format is used in the assembly of nonlinear system models. The
Volterra dynamic representation was recognized and defined for first time in this work.
The Volterra dynamic representation is important in the MMM because this external
port based model format protects internal design details. This procedure generates
the Volterra dynamic port-based system models in an explicit form for the first time.
The method to obtain frequency domain port-based Volterra system models from
MIMO port-based ODEs is discussed in Chapter 4.

The method to assemble physical MIMO external port—based Volterra system
models is an original contribution to mechanical engineering because it is a direct
procedure to assemble nonlinear, port-based, Volterra system models using port-
based physical constraints. Volterra representations are appropriate to the MMM
approach to distributed model assembly because they protect internal design details.
This method extends the MMM approach to nonlinear port—based models through a
method that is computationally efficient because it is recursive and does not require
matrix inverses. These properties make this Volterra-based assembly procedure ideal
for distribution and assembly of nonlinear physical engineering system models. The
method to assemble physical MIMO port-based Volterra representations is described

in the chapter 5.

1.3 Limitations of the Nonlinear MMM Extension

The nonlinear Modular Model distribution and assembly methodology has two (2)

limitations. These limitations are:

1. The nonlinear subsystems must be modeled using port-based Volterra representa-

tions.

2. The assembled system model is valid if the subsystem models used in the assembly
are valid.

The nonlinear subsystems must be modeled using port—based Volterra representa-
tions because the formats used in the nonlinear MMM extension developed here are
port-based Volterra models. Not all nonlinearities can be represented using Volterra
models. A non-linearity can be represented using a Volterra model if it is analytic
or if it has a fading memory. The analytic requirement for Volterra models was first
recognized by Brilliant [11] and later by Sandberg [12]. Because many common non-
linear systems are modeled with non-analytic non-linearities, Boyd [13] reduced this
limitation by showing that for non—analytic nonlinear operators with fading memory,
a Volterra representation can be obtained for a useful set of bounded input signals.
Although port-based Volterra models are not possible for all nonlinear systems, SISO
Volterra models are common in engineering literature and have been used for a variety
of engineering applications [14],[15] and [16].

The assembled system model is valid if the subsystem models used in the assembly
are valid. In the MMM process it is assumed that each provided Volterra subsystem
model is valid at a specific input-output operating condition. When these subsystem
models are assembled using the MMM procedure, the input-output operating condi-
tions for each subsystem in the assembly do not change, then the resultant assembled
model is valid. Component operating conditions do not change because the MMM
method provides an explicit method for computing the operating point inputs of the
assembly required to Operate the components at the desired input-output conditions.
This explicit method is explained in the Chapter 3 of this thesis. It is important
to clarify that the validity of the Volterra subsystem models is not address in this
work. This work only deal with the problem of assembly and distribution of nonlinear

models.

CHAPTER 2

Background

2.1 Internet Distribution of Models

The Internet has enabled a global engineering design strategy where collaborative de-
sign teams perform irrespective of physical distance. Models of components provided
from different sources must be distributed and assembled to generate system models.
This section will discuss the functionality of the networked software agents required
to implement the model information flow.

The performance of a networked model distribution system (Figure 2.1) relies on
the operation of autonomous and flexible computational systems called Agents [17].
Four agent classes are used: Component agents, Assembly agents, the Agent Reg-
istry and the Query Ontology. Initially an User connected to the Internet, consults
to an Agent Registry for the locations of Assembly and Component Agents on the net-
work. The User uses standardized queries whose format is obtained from the Query
Ontology which publishes an organized list of these queries. An assembly agent assem-
bles models by using queries to lower level components. Assembly and Component
agents, respond to queries with models in the standard model format specified by the
Ontology.

The networked distribution of models is hierarchical. This characteristic is the

result of the agent’s capabilities to perform either as a client, as a server or both.

 

 

   
   

 

 

 

 

Component 1 ] Query
Agent Ontology

 

 

 

 

l

 

 

 

Com onent 2
I/ Assembly ‘~\ X
’1’ \\ gent
”’ I \ “ \
-_4{.’ _____ ..1.’.I.__~ f \\> \E._

 

 

( f ‘
E Component 1 mt; Component 2 E
A) ‘ as I_

A

Figure 2.1. Physical View of a Networked Modeling System

 

 

 

 

 

 

 

 

 

 

 

 

Server

Component Agent

 

tier 2 Component Agent

I Server I I Server I
tier 3

 

 

 

 

 

 

 

 

Component Agent Component Agent Component Agent

 

 

 

Figure 2.2. Logical Information Flow View of a Networked Engineering System.

The User can perform only as a client, requesting model information from lower-level
agents. Assembly agents can perform either as clients requesting model informa-
tion from lower-level agents or as servers providing model information to higher-level
agents. Finally, component agents can perform only as servers, providing model in-
formation to higher-level agents.

In a three level network example (Figure 2.2), the User requests model informa-
tion from the server in the tier 1 assembly agent. The tier 1 client requests model
information from servers in tier 2 agents. The tier 2 client requests models from tier 3
component agent servers. Starting at the lower tier in the system, model information
is provided by servers to higher-level assembly agents that assemble them into yet
higher-level assembly models. These assembled models are then provided to agent

servers that respond to queries from above.

2.2 The Standard Linear Physical Model Formats

A physical system is an entity separated from the environment that interchanges
energy through a boundary [7]. Physical systems are composed of interacting compo-
nents that perform in a synchronized way to generate an energy flow. This energy
flow is transferred through physical connections consisting of input-output pairs called
ports [18]. The product between the input and the output variables of a port defines
the energy flow through the port. Positive energy flow through a port is defined as
the work done on the system. The total energy flow in a physical system is the sum
of all the energy flows through each of its ports.

Physical systems can be modeled using a port-based approach. Port-based models
always have an equal number of inputs and outputs because an input-output pair
defines each port [7]. Port-based models are considered external models if they are
given as functions of external port variables. A valid external port-based model has

independent external work ports. This characteristic requires the number of model

11

equations to be equal to the number of system ports outputs. An external, time-
invariant, port-based dynamic linear physical system model having r external ports

has r equations in the form

N0y(t) +---+N,- ((11:29)) +---] — [M0y(t) +---+M,- (£1330) +--] (2.1)

 

 

where, N,- V( 0 S 2' S n) is a (r x r) matrix of constants, Mj V( 0 S j S m) is a (r x r)
matrix of constants, y(t) is the (r x 1) system output vector and u(t) is the (r x 1)
system input vector. Causal physical system models require m _<_ n [19]. Applying

the Laplace Transform to the linear ODE (2.1) and factoring common terms,
[N050 + - - - + ann] Y(s) = [M030 + - - - + Mmsm] U(s) (2.2)

where Y(s) and U(s) are respectively the Laplace transform of the physical system
output y(t) and input u(t) . Defining the two polynomial matrices

N(s) —_— [N030 + le + - - - + ann]

(2.3)
M(s) = [M030 + M15 + - - - + Mmsm]
and substituting them into (2.2), yields
N(8)Y(S) = M(8)U(8) (24)

where N(s) and M(s) are (r x r) matrices of polynomials. Two model representations
can be derived from (2.4). These are the dynamic model representation [20] and the
transfer function model representation [21].

The dynamic model representation
P(s)Y(s) = U(s) (2.5)

uses the (r x r) matrix

P(s) = [M(s)]‘1N<s) (2.6)

This first model representation is a dynamic extension of the traditional format used

in dynamic structural analysis [5]. The dynamic matrix P(s) exists only if M(s) is

12

nonsingular. When M(s) is nonsingular, the effects of the port inputs in the equations
are independent. Model representation (2.5) exists even if the matrix N(s) is singular,
allowing models where the effects of their outputs in the equations are dependent.

The transfer function model representation
Y(s) = G(s)U(s) (2.7)

uses the (r x r) matrix

G(8) = [N(5)l—1 M(S) (2-8)

This second representation is the traditional format used in system dynamic analysis
[22]. The transfer function matrix G(s) exists only if N(s) is nonsingular. When N(s)
is nonsingular, the effects of the outputs in the equations are independent. Model
representation (2.7) exists even if the matrix M(s) is singular, allowing models where
the effects of their inputs in the equations are dependent.

The linear state space model representation [22],

x=Ax+Bu

y = Cx + Du (29)

is another commonly form used to obtain the external model representation (2.4). A
state space model representing a physical system with r ports and p internal states
uses a (r x 1) output vector y(t), a (r x 1) input vector u(t) and a (p x 1) state vector
x(t). Space state models are internal representations because they are functions of
internal state variables. The internal state variables can be removed to obtain an
external representation in the form (2.4). The transfer function model representation

can be derived from (2.9) by using the well-known equation
Y(s) = [c [51 — A]“1 B + D] U(s) (2.10)

The equivalent model format (2.4), using equation (2.10) is
I|[sI —— A]| = [C adj(sI —- A) B + [[31 — A]|D] U(s) (2.11)

13

where adj is the matrix adjoint operator. Comparing (2.4) and (2.11),
N(s) = [[31 — A“ (2.12)

M(s) = [Cadj [sI — A] B + ”81 — A]|D] (2.13)

In general, external model representations derived from (2.1) can describe systems
with outputs that produce dependent effects on the port equations (the matrix N(s)
can be singular). In contrast, the external model representations derived from the
state space representation (2.9) do not allow this characteristic because the resultant
matrix N(s) is always nonsingular since N(s) = [51 — AI 75 0 .

Port-based physical systems can be modeled by the dynamic model representation
(2.5) or by the transfer function model representation (2.7 ) If the matrices M(s) and
N(s) are nonsingular, both model representations can be obtained either by using
the external model (2.4) or by using the state space model (2.9). These two model
representations are the standard formats used by the linear MMM for distribution
and assembly of models.

In the next section, the two properties that characterized the assembly of physical
systems are introduced. These properties are represented mathematically by two
constraints equations. These equations are used in the linear and nonlinear MMM

assembly method.

2.3 Physical Assembly Constraints

Physical model assembly is characterized by two conditions. The first condition is
associated with the word assembly. The assembly of two or more system models means
their response or outputs variables are equal at the geometric location of the physical
assembly. This condition will define the relationship between port output variables.
The second condition is associated with the word ”physical”. All physical system

connections must satisfy conservation of energy. The assembly of physical systems

14

requires both conditions at physical connections: equal system output response and
energy conservation.

The connection of physical subsystem models is executed by joining port variables
that exchange energy. The output and input pair of each port is determined using the
two concepts characterizing the assembly of physical systems. The output variable is
selected as the variable physically constrained to be equal to other output variables
when ports are connected. The input variable is selected, as the port variable required
to compute units of energy when it is multiplied by that port’s output. Energy

conservation requires all connected ports to be in the same energy domain.

 

r Energy Domain | Output (y) ] Input (u)J

 

lOIl OI'CB
1011

O ume
O ume
11X

 

Table 2.1. Input-Output MMM Causality for Different Energy Domains

Input and output port variable pairs for different energy domains are shown in
Table 2.1. Connected port outputs measured with respect to the assembly’s system
reference are equal. Ports assembled in the electric energy domain require equal
potentials. Ports assembled in the mechanical energy domain require either equal
angular or linear displacement. Ports assembled in the hydraulic or acoustic energy
domain require equal pressure. Ports assembled in the heat transfer energy domain
require equal temperatures. All these required conditions describe physical assembly
in their respective energy domain. The associated port inputs in each energy domain
permit the computation of energy. This standard facilitates a the model assembly
process. The naturally constrained port variables allow the use of standard physical
assembly constraints to assemble models.

Joins are used to enforce these assembly constraints. A join is not a model of

15

a physical connection subsystem and does not store nor dissipate energy. Joins are
mathematical mechanisms that provide the proper physical connection constraints
for connecting the ports of physical components within a single energy domain [8]. A
join is graphically represented as a circle enclosing the letter J for Join (Figure 2.3).
The lines represent ports with the direction of positive energy flow indicated by the
arrowheads.

Two kinds of control volumes are defined in any assembly. The component control
volume (c) interchanges energy with the assembly control volume (a) through the
component output yc(t) and input uc(t) port variables. The assembly control volume
interchanges energy with the environment thorough the assembly output ya(t) and

input ua (t) port variables.

 

\
\
‘\

 

\
“
‘\

 

 

 

 

—
”- -‘

 

 

 

: -. ‘
Physical ] |
System ,u' ’

I J I, i”:

 

 

'-_---_
’0

I,

A
H
__ 4’
ID
tuetuuonAug

\\ “‘\ I; I I,
\ Physical //
\\\\\‘\ System ,x/ ” (ya (0: “a (1))
\\‘;: ___________ “I” ,x” ( (t) (t)
\\\“ :::: /’,’ y c , “c )

a
’
"
-————-'

Figure 2.3. Unconstrained and Assembly Control Volumes

A simple relationship determines the number of assembly variables. If a set of
components with r port variables uses f connections to connect p component ports,
the assembled system has l = r+ f — p pairs of port variables. In a practical assembly
f _<_ p and l S r . The dimension of the component and assembly port variable
vectors are respectively (r x 1) and (l x 1).

The two conditions that characterize a physical systems assembly can be repre-

16

sented using two equations. The first equation constrains connected outputs to be

equal. For a join connecting d outputs,

3110‘) = --- = u(t) = -°-= yd(t) = ya(t) (2-14)

where yz-(t) is the ith connected component output at the join and ya(t) is the output
assembly variable at the join. For both joins in (Figure 2.3), d = 2. For a set of l
joins generating l assembly output variables from r component ports, the relationship

(2.14) between the r component outputs and the l assembly outputs is,

yc(t) = Sya(t) (2-15)

where S is a (r x l) matrix of constants called the constraint matrix. The second
equation, constrains the energy stored at all the connections to be zero when the
work done internally on the component ports equals the external work done on the

assembly,

ufit>yc<t> = uiitiyam (2.16)
This equation states that the total energy supplied to the component port equals
the total energy supplied to the assembly port. The equation (2.16) includes input

and output variables and is difficult to use. An equation easier to use relates only

component and assembly inputs and can be derived substituting (2.15) into (2.16),
uZit>Sya<n = uZ<t>ya<t> (2.17)
From equation (2.17) it follows that,
STuc(t) -.-. ua(t) (2.18)
Applying the Laplace Transform to (2.15) and (2.18),
Yc(s) = SYa(s) (2.19a)
STUC(s) = Ua(s) (2.1%)

17

Equations (2.19a) and (2.1%) are the two constraints used to assemble physical mod-
els. This two constraints are satisfied by any assembly, independently if the compo-
nents are linear or nonlinear.

The assembly of three components is shown (Figure 2.4). These three components
have a total of eight (r = 8) ports. Component 1 has three ports: 1,2 and 3. Com-
ponent 2 has 1 port: 4. Component 3 has four ports: 5,6,7 and 8. Five additional
ports (a, b, c, d and e) are defined for each of the five (f = 5) joins of the assembly
(J1, J2, J3, J4, and J5). These joins are used to connect eight (1) = 8) component
ports. The total number of ports for the assembled system is five: (I = r + f — p = 5).

They are the ports a, b, c, d and e.

e

Y2=ye ‘ ys=ye
yfi yc

. 8 1‘ 5 ’
- .4>»5:“ ' . y7=yd

‘ y4= yb
Y3 =yb ys =yb

b

 

 

 

 

any.

 

In.

 

at. ' “a.
1H"
”I

 

 

 

 

 

 

 

 

 

 

Figure 2.4. Assembly of Three Components Through Five Connections

A constraint is defined for each of the eight connected component port outputs.
At the join (J1), the output of component port 1 is constrained to be equal to the
output of assembly port a. At the join (J2), the outputs of component ports 3, 4

and 5 are constrained to be equal to the output of assembly port b. At the join (J3),

the output of component port 6 is constrained to be equal to the output of assembly
port c. At the join (J4), the output of component port 7 is constrained to be equal
to the output of assembly port (I. At the join (J5), the outputs of component ports
8 and 2 are constrained to be equal to the output of assembly port e. The matrix

form of the equations that relate component outputs and assembly outputs is,

 

 

 

 

 

 

’y11'10000'

y2 00001~-

3,3 01000 3;:

:: =31333 ,
ya 00100[yd

3,7 00010 ye,

[3,8, L00001_

This equation is an example of (2.15) that defines the constraints between component
and assembly outputs through the constraint matrix S. A constraint equation is
defined for each of the five joins used. At every join the sum of the component port
inputs must equal the input of the respective assembly port. For joins (J1), (J2) and
(J3) the result is trivial because only a single component port is joined. At the join
(J2) the sum of the inputs of component ports 3,4 and 5 must equal the input of the
assembly port b. At the join (J5) the sum of the inputs of component ports 2 and
8 must equal the input of the assembly port e. The equation that relates component

port inputs and assembly port inputs is,

 

 

 

 

.1“.
'10000000'“2 'ua‘
00111000 ”3 ub
0 0 0 0 010 0 “4 = uc (2.20b)
00000010 “5 ud
00000001 “6 ue
- -W _-
-118.

 

 

The matrix above equals the transpose of the constraint matrix S in (2.20a). The next

section shows how the two physical dynamic constraints equations (2.19a) and (2.19b)

19

are used to assemble external physical system models. Additionally, it is shown that
two different model representations, one model representation used for distribution
and assembly and the other model representation used for simulation and analysis

are required to effectively implement the MMM process.

2.4 The Modular Modeling Assembly Procedure

The Modular Modeling Method (MMM) uses a systematic process to assemble dy-
namic matrix models (2.5). These standard format models have port pairs standard-
ized through the two concepts that generate constraints (2.19a) and (2.19b). It is
shown in this section that the dynamic matrix representation is a convenient form
for assembling physical system models from component models. The MMM assembly
includes four steps. These steps are: 1) Generation of the unconstrained component
model, 2) Generation of the constraint equations, 3) Generation of the assembled
model and 4) Generation of the condensed model.

The generation of the unconstrained component model is the first step. This
process formulates a matrix of component dynamic models (2.5) in diagonal form.

An assembly of k component models with a total number of r ports yields,

P1(S) 0 0
Pc(s) = 0 P209) 0 (2.21)
_ 0 0 O Pk(s) .

 

 

where Pi(s), (1 S 2' _<_ k) is the (r,- x r,) dynamic matrix of the ith component and
Pc(s) is the (r x r) unconstrained matrix constituted by the dynamic matrices of all

the assembly’s components. The number of ports in an assembly of k components is,

k
r = Zri (2.22)
i=1

The unconstrained component model relates the (r x 1) component output vector YC

2O

to the (r x 1) component input vector UC through the equation,
Pc(s)YC(s) = Uc(s) (2.23)

The generation of dynamic constraint equations in format (2.19a) and (2.19b)
is the second step. This process builds the constraint matrix S that relates uncon-
strained component variables to constrained assembly variables. The matrix S is
dependent on the connection topology between components. This matrix equates
each component output in Y0 to an assembly output in Ya based on the constraints
associated with the assembly. An example is shown in the next section.

The generation of the assembled model is the third step. This process is executed
by applying output (2.19a) and input (2.19b) constraint equations to (2.23). Initially,

(2.19a) is substituted into (2.23) to yield,
Pc(S)SYa(8) = Uc(8) (224)
Multiplying both sides of (2.24) by ST and using (2.19b) yields the assembled model,
name) = U.(s> (2.25)
where
Pan) = STPC(s)S

is the (l x l) assembly dynamic matrix and Ua(s) and Ya(s) are the (l x 1) assembly
input and output vectors.

The dynamic model representation (2.5) is particularly effective as the standard
assembly MMM format because constraints (2.19a) and (2.19b) can be easily applied.
Additionally, the overall process is recursive since the assembled model (2.25) is again
in the standard format (2.5). In contrast, the transfer function model (2.7) cannot be
used as the assembly MMM format because (2.19a) and (2.19b) cannot be substituted

directly. This process requires the (r X I) left inverse of the matrix S. This matrix

21

does not exists because in the matrix S the number of rows is greater than the number
of columns. Even though standard format (2.5) is effective for the model assembly
process, it is not appropriate for simulations because is not possible to use it to solve
directly for the outputs given the inputs. A different model format is required for
simulation.

The transfer function model representation (2.7) is chosen as the MMM simulation
format because it can directly be used to solve for the outputs given the inputs. This

format is obtained by inverting (2.25) to yield,
Ya(8) = Ga(5)Ua(8) (2-26)
where
Gas) = [P.(sn‘1

is the transfer function matrix. Model (2.26) is the traditional model format used in
the analysis and simulation of dynamic systems [22].

The two different MMM system representations (2.25) and (2.26) are important.
The analysis above demonstrates that the MMM modeling format (2.25) is partic-
ularly well suited to recursive assembly of models using output constraints but not
appropriate for perform model simulation. The analysis above also reveals that the
model format (2.26) is particularly well suited for simulation but not appropriate to
perform recursive system model assembly. If only one format is used, the advantages
of both model representations are not obtained. MMM used the model representation
(2.25) for distribution and assembly and the model representation (2.26) for analysis
and simulation. This is a clear difference with the traditional methods that use the
transfer function model format (2.26) as the unique model representation.

The assembly model (2.25) is formulated in terms of assembly port variables con—
sisting of both internal and external assembly ports. Internal ports are those ports

whose inputs are zero during normal use of the model. These zero input ports can

22

be removed from the model. Additionally, these ports might be removed to protect
the proprietary response of the internal model ports. Conversely, external ports must
be accessible to the users and should remain in the model. A process to eliminate
the internal ports is required. This process is called condensation and the resultant
model is only a function of external port variables.

Model condensation is the fourth and final step. This process uses a symmetric
transformation matrix T to reorganize assembly variables into internal and external

varaible vectors in the form,

Ye(s) __
T ] “(8) ] _ Ya(s) (2.27a)
Ue(8) _ T 5
[ U.(s) ] _r Ua( ) (2.27b)

where the port input Ue(s) and output Ye(s) are external and the port input U,(s)
and output Y,(s) are internal. To remove internal port variables substitute (2.27a)

into (2.27b),

 

Ye(s) _
Pa(s)T Yi(s) ] — Ua(s) (2.28)
Multiplying both sides of (2.28) by TT and substituting (2.27b),
T Ye(5) __ Ue(5)
T P“(S)TlY.-<s> l ‘ l Ui(5)l (2'29)
Defining the matrix,
T _ Pee(8) Pei(3)
T Pa(S)T _ i Pie(3) Pii(8) i (230)
yields the form,
1368(5) Pei(3) ] [Ye(3) ] = [ [16(3) ]
lees) p,,(.) Y.(s> U.<s) (2'31)

where Pee (s), P6,;(s) , P,e(s) and Piz-(s) are dynamic matrices. Because internal
assembly inputs will not be accessible in the condensed model, these inputs are set to
zero (Ui(s) = 0) to solve for the external inputs as a function of the external outputs

in the form,

Pe(s)Ye(s) = Ue(s) (2.32)

23

where,

Pe(5) = P€€(S) _ Pei(3)Pii(5)P;.-31(3) (2-33)

is the condensed dynamic matrix. The condensation process is recursive because the
condensed model is again in the standard format (2.5). Internal dynamics matrix

inverse P31 is always computable (See Appendix A).

2.5 Example of a Linear MMM Assembly

The assembly of a mechanic transmission linear model and an electric generator linear
model is developed in this section. The first component is an electric generator with
two ports (Figure 2.5). The generator’s first port variables are the input electrical
charge (11,1 (t) and the output voltage potential e1,1 (t) . The generator’s second port
variables are the input rotational torque r1,2 (t) and the output angular displacement
612 (t). The second component (Figure 2.5) is a mechanical transmission with two me-
chanical ports. The transmission’s first port variables are the input rotational torque
r2,1(t) and the output angular displacement 62,1(t) while the transmission’s second
port variables are input rotational torque 72,2 (t) and the output angular displacement
622 (t) In total, the two components have r = 4 ports.

The assembly is executed thorough three (f = 3) joins that connect three (l = 3)
assembly ports to the four (r = 4) component ports. The left joint connects the ex-
ternal electric charge-output potential port pair (ea,1(t), qa’1(t)) to the internal port
pair (€1,1(t), q1,1(t)) . The right joint connects the external torque-angular displace-
ment port pair (60,3(t), ra,3(t)) to the internal port pair (62,2(t), r2,2(t)) . Finally the
middle joint connects the internal port pairs (61,2(t),7'1,2(t)) and (62,1(t), 72,1(t)) to
the external torque—angular displacement port pair (0a,2(t), Ta,2(t)).

The dynamic model representation for the generator is

P1,11(8) P1,12(8) E1,1(s) _ Q1,1(8)
P1,21(5) P1,22(8) l [ 91,2(3) ] _ [ 112(3) ] (2-34)

24

(10,10)

 

]ea,1(f) 10,2(1) ”9.130) Ta,3(t)

 

1490,30)

 

 
    
 

‘— €1,1(t) 61,2 (9+ 902,10)

 
   

q, (t)
_+ 11 mo;— 72,10, 922i),

72,2(1)

 

 

 

Figure 2.5. Assembly of a Transmission and an Electric Generator

where E1,1(s), 6214(3), 913(3) and T1,2(s) are respectively the Laplace trans-
form of the time variables €1,1(t), q1,1(t), 013(3) and 71,2(t), the complex variable
functions P1711(s), P1,12(s), P1,21(s) and P1,22(s) are ratio of polynomials respec-
tively representing the relations between the input-output pairs (Q1,1(s),E1,1(s)),
(713.2(3),E1,1(S)),(Q1,1(S),912(3)) and (T1,2(8),91,2(3))-
The transmission dynamic model representation for the is
P2,11(3) P2,12(8) ] [ e2,1(8) ] = [ T2,1(3) ] (235)
P2,21(8) P2,22(5) 92,2(8) 722(5)
where 82,1(5), T2,1(s), 82,2(3) and T2303) are respectively the Laplace trans-
form of the time variables 92,105), 7271(t), 02,2(3) and T226), the complex
variable functions P2’11(s), P2,12(s), P2121(s) and P2,22(s) are ratio of poly—
nomials respectively representing the relations between the input-output pairs
(721(8),92,1(8)),(T2,2(S),92,1(8)):(T2,1(8),922(5)) and (722(3): 922(3))-
The construction of the unconstrained component model is the first step. The
unconstrained matrix is built by inserting in the main diagonal of a square matrix,

the dynamic matrices of component models (2.34) and (2.35). The other elements of

25

the unconstrained matrix are set to zero.

136(8)

"P1,11(8) P1,12(8)
P1,21(8) P1,22(8)
0 0
0 0

 

_

0
0

P2,11(8) P2,12(8)
P2,21(S) P2,22(5).

0
0

The unconstrained component model Pc(s)Yc(s) = Uc(s) is,

FP1,11(8) P1,12(8)
P1,21(8) P1,22(5)
0 0
0 0

 

0 0

0 0
P2,11(8) P2,12(8)
P2,21(3) P2,22(S)_

 

' E1,1(8) I
91,2(8)
92,1(8)

. 92,2(8) .

 

 

 

P Q1,1(8) '
T1,2(3)
T2,1(3)

 

 

t T2,2(3) 1

(2.36)

(2.37)

The generation of the constraint equations is the second step. Initially the time

variable equations that relate component output to assembly outputs are written.

 

 

 

P (31,1(t)1 I.1 O 0 1
(91,2(15) _ 0 1 0 30,18;
92,1(t) — 0 1 0 gait
_92203), _0 0 1_ “’30
1 0 0‘T q 10) 'q1,1(t)‘
0 I 0 Ta,2(t) : T1303)
0 1 0 T“, (t) mu)
-0 0 ll 0’3 .7220).
Applying the Laplace Transform to (2.38a) and (2.38a),
E1,1(S) 1 0 O
ems) 0 1 0 gal]:
92,1(3) 0 1 0 9“,?
923(3) 0 0 1 03(3)
'1 o 0"" Q 1(3) "621.1(5)
0 I 0 Ta,2(8) : 712(3)
0 I 0 Ta, (8) 7121(8)
_0 0 1_ “’3 run)
The constraint matrix is _ q
1 0 0
0 1 0
S— 0 1 0
_ 0 0 1_

 

 

 

 

 

 

 

 

 

 

 

(2.38a)

(2.38b)

(2.39a)

(2.3%)

(2.40)

The assembly model generation is the third step. Using (2.39a) and (2.39b) to gener-

ate the assembly’s dynamic matrix Pa(s) = STPCS,

P1,11(8) P1,12(8) 0 Ea,1(8) Qa,1
P1,21(8) P1,22(8)+P2,11(8) P2,12(3) e(1,2(8) = Ta,2 (2-41)
0 P2,21(3) P2,22(8) 9a,3(3) Ta,3

Equation (2.41) is the dynamic model of the generator-transmission assembly.
This model in the standard format (2.5). At this point, any port of (2.41) can be
selected as internal or external. Internal ports make internal variables unavailable to
the user. Because internal ports are not accessible, internal inputs are set to zero.
The condensed model is determined assuming that the internal variables are 90,2(3)
and Ta,2(s) and the external variables are Ea,1(s), Qa,1(s), 603(3) and Ta,3(s).

The condensed model generation is the fourth step and uses the matrix,

100
T=001 (2.42)
010

to reorganize the assembly variables into external and internal vector variables,

Ea,1(5) Ea,l
T ea,3(8) = 902 (2.4321)
e(1,2(3) ea,3
Q1,1(S) Q1,1
Ta,3(s) = TT Tag (2.43b)
Ta,2(5) Ta,3
Using the definition in (2.30), (2.31), (2.32) and (2.33) yields the condensed model,
Pee(3)P ei(3) ] [in (1 ,1(5) ] : [Q1,1(3) ] 2 4
l as) Pas) 3e ) 123(8) ( ' 4)

where

P1.12(S)P1,21(8)
Pee(8) = P1.11(5) ’ P1,22(3)+P2,11(3l

_ P1,12(8)P2,12(8)
Pei“) _ _P1,(22 (§l,+P211(8)

 

 

 

. (8———) P2,21(3)P1,21(3)
“3 P1, 22(8)+P2, 11(Sl

P2,21(3)P2,12(S)
1022(5) 2 P2722“) — P1,22(S)+P2,11(8)

 

27

The condensed model is also in the standard format (2.5) and can be assembled to
other standard models using the same algorithm. The condensation process is again
recursive.

The analysis and simulation of the generator-transmission system assembly, re-
quires the simulation format (2.26). This model is obtained by inverting the con-
densed model (2.44) to yield,

[ Ea,1(5) ] _ ] Gee(8) Geifs) ] [ Q1,1(8) ]

903(5) — Gte(8) (322(3) Ta,3(s) (2'45)

Where,
Gee(8) = [P2,22(S)P1,22(S) + P2,22(3)P2,11(8) — P2,21(8)P2,12(8)] / 17(8)
062(3) = [P1,12(8)P2,12(S)l / 17(3)
0212(8) = [P2,21(8)P1,21(3)l / 13(3)
022(8) = [(P1,22(8) + 1211(8)) P1,11(S) — P1,12(8)P1,21(8)l / 17(8)
0(8) = Gee(8)P1,11 - P1,12P1,21P2,22

2.6 Summary

In this chapter, the recognition that two model formats are necessary for physical
system model assembly and simulation is addressed for the first time. The use of two
model formats decreases the computational time required to assemble and to simulate
physical system models. These two model formats are the transfer function represen-
tation and the dynamic representation. The transfer function representation is the
traditional format used for simulation. This thesis shows, that in a constraint-based
model assembly, the Transfer Function representation is not appropriate for MMM
assembly because a matrix inversion is required. Matrix inversion demands substan—
tial computational time, a clear disadvantage in a global distributed environment.
Additionally, for this inverse to exist, system models are limited to have independent

inputs and independent outputs. This is a big limitation on allowable models that

28

excludes many dynamic system models. In contrast, if the format used for assembly is
the dynamic representation, no matrix inversion is required and system models with
dependent inputs and dependent outputs can also be assembled.

This chapter also introduces the process of model condensation. Model condensa-
tion is the removal of ports that are considered interior to the model. Removing those
degrees of freedom reduces the overall size of the model, and additionally, protects the
proprietary information of the model. Condensation of models significantly increases
the difficulty associated with reverse engineering a model. Invariably the number of
internal model parameters exceeds the number of inputs/ output pairs making reverse
engineering difficult. In the next chapters the MMM process for assembly physical

nonlinear models is presented.

29

CHAPTER 3

Assembly of Affine Physical System
Models

3.1 Affine Physical System Models

Affine systems are important in the MMM procedure because the assembly of general
differentiable nonlinear physical models performing at a constant operating point can
be executed as the assembly of its affine approximations. In an affine system the
inputs and outputs exhibit a proportional relationship, but the model outputs are
nonzero at zero input. Affine systems are nonlinear because they do not obey the
two linear system properties: superposition and homogeneity [23]. Affine systems
can be static or dynamic. A static affine system is characterized by an input-output
relationship that is independent of the system’s input and output time derivatives.
Static affine systems are discussed first, then, dynamic affine systems are presented.

A static affine system having p outputs and q inputs is,
y = Ku + c (3.1)

where y is a (p x 1) vector that represents the system’s outputs, u is a (q x 1) vector
that represents the system’s inputs, K is a (p x q) matrix of constants and c is a
(p x 1) vector of constants called‘the bias vector. The static affine model (3.1) is

non-linear because it does not obey superposition and homogeneity. Superposition is

30

not satisfied by (3.1) because given two arbitrary inputs 111 and Hz the relation:

Y(u1 + u2) = Y(u1) + Y(u2) (3-2)

does not hold. Homogeneity is not satisfied by (3.1) because given an arbitrary

constant A and any input 11, the relation

y(Au) = M(U) (33)
does not hold.

Two examples of physical affine systems are an electrical battery model with non—
zero open circuit potential (Figure 3.1) and a mechanical system with fixed load
(Figure 3.2). Both systems have input-output models with non-zero output at zero

input. The static affine electric battery model (Figure 3.1)
e = —R2' + V (3.4)

has output potential 6, input load current i and battery internal resistance R. Here

the bias constant V, is the battery potential at zero current.

 

 

 

 

 

 

A

 

 

~ Current

Figure 3.1. Electric Battery

The dynamic affine pulley-load assembly model (Figure 3.2) is described by the

nonlinear ordinary differential equation

167+ k6 = 7' + WR (3.5)

31

with output angular displacement 0, input torque 7, spring parameter k , mass mo-
ment of inertia J, and bias constant WR. The bias constant depends on both the

pulley radius R and the load W. At static equilibrium model (3.5) becomes,

where, Oeq is the angular displacement at equilibrium for a given constant torque r.
The static model (3.6) does not include the point (06g, r) = (0,0) . Like the previous
example, the static (3.6) and the dynamic (3.5) models are nonlinear because they

do not obey superposition and homogeneity.

 

7’
I
I
I
I
I
5
I
y
I
I
I
I
I
I_

 

 

Figure 3.2. Pulley-Load Assembly

An affine system is an intermediate result of any linearization process. A lineariza-
tion (Figure. 3.3), usually generates a linear model with respect to the operational
variables but an affine model with respect to the real physical variables. Physical
variables are most appropriate for use in a model’s assembly process because when
models are connected to form a system model, their behaviors are constrained [6] and
these constraints are in function of physical variables.

An affine approximation of the external nonlinear model,

f(vaiyiu°iu7fiifii"°):0 (3.7)

32

 

 

 

 

 

 

Figure 3.3. Linearization Yielding and Affine System

where f(o) is a (r x 1) vector of functions, y(t) is the (p x 1) system physical output
vector and u(t) is the (q x 1) system physical input vector, can be obtained through
a Taylor expansion in five steps.

In the first step, all the input and output derivatives of (3.7) are set to zero,
f"'(y,u)=f(y,0,---,u,0,-~)=0 (3.8)

where for any y on the map, at least one u exists.

In the second step a point (y, u) on the map (3.8) is selected. Traditionally, a
desired value of the output y 2 yo is defined first to determine one input 11 = uO
that satisfies,

f*(y0, u0) = 0 (3.9)
In the third step, the equilibrium point

"zuo’f°:i’0=m=0 (3.10)
y=y0, y0=y0=---=0
is selected as the system’s operating point (op). Applying Taylor series expansion to
equation (3.7) around the operating point (3.10)

f(y, 11. - - -) = f(y0. no, - - -) + Dy(f)l0,,(y - yo) + 211D§(0|0p(y — yo)2

3.11
+ . . . +Du(f)|0,,(u- “0) + %,Dg(t)|op(u_ u0)2 +... = 0 ( )

33

where each element of the matrix D§(h(g, g)), - - - )) is the the kth partial derivative
of the vector function h respect the vector g. Since (3.10) is a solution of (3.7) the

first right side term of (3.11) is,
f(y0,u0,0,---) :0 (3.12)

In the fourth step, equation (3.10) and (3.12) are substituted into (3.11) and

retaining only the 1“ partial derivative terms,
f(y, u, . .-) a“ Dy(f)|0p(y - yo) + D5,(f)|0py+ Dymlopf' + . ..
+Du(0lop(u — “0) + D,-,(t)|opu + Dfimlopa + . ..
Equation (3.13) can be reorganized in the form
f(y,u, . .. ,y, u) 9.: Dy(t)|opy + Dy(t)|0py + D5,(f)]0py+ . ..
+Du(0i.,u + Du(fl|0,,fi + Damiopu + - ~ (3.14)
-Dy(0lopy0 - Duffllopuo

All the partial derivatives evaluated at the operating point (op) in (3.14), yield con-

(3.13)

stants matrices. The fifth and final step is to define:

1. n (r x p) matrices of constants (1 gig n) Ni: [Bf/6(diy/dti)]|0p

2. m (r x q) matrices of constants (1323 m) M,- = [3f/8(diu/dti)]|0p

where n and m are respectively the highest derivative of the output and the input of
equation (3.7).

3. A (r x 1) vector of bias constants c=- [Bf/8y] [opy0 — [Bf/au] Iopu0

to write (3.14) in the form,

f(y,u,--- ,y,u)§N0y+N,y+N2y+...+

. .. (3.15)
M0u+M,U+M,U+---+c

Equation (3.15) evaluated in its physical variables is affine if the vector c does not

vanish. The vector c will vanish if,
[Bf/0y] [01,,y0 + [Bf/Bu] lopuO = O (3.16)

Operating points that satisfy condition (3.16) can be approximated in its physical

coordinates by a linear model and are a special case. This fact would imply either

34

that the operating point Of the assembly is the origin (yo, no) = (O, 0), what in most
cases is not practical, or that even though the Operating point is not the origin, the
linearized system intersects the origin. Any differentiable nonlinear system that does
not satisfy condition (3.16) can be locally approximated in its physical variables by
an affine model (3.15).

An affine model describing a physical system with r external ports, requires r

equations in the form,
N0y+N1y+N2y+---+M0u+Mlu+M2ii+---+c=0 (3.17)

where respectively V(O S t S n) and V(0 g j S m), N, and M j are (r x r) matrices
Of constants and y(t), u(t) and c are (r X 1) vectors. Physical systems are causal

requiring [19]. The bias vector,

c = —N0yO — MOuO (3.18)
substituted into (3.17) yields,
N0y+N1y+N2$r+m+M0u+M1u+M2ii+---—N0y0 —M0u0 = 0 (3.19)

A procedure to assemble affine physical models in their physical variables is proposed
in this chapter. The procedure satisfies the four characteristics required by a global
engineering strategy and applies to any differentiable nonlinear physical model. The
chapter is organized as follows: Affine systems were defined in the first section. The
networked distribution process for affine models is explained in the second section.
The standard model format is described in the third section. The constraints required
to assemble affine physical models are defined in the fourth section. The MMM
algorithm to assemble affine physical models is shown in the fifth section. An example

is shown in the sixth section. Finally conclusions are provided.

35

3.2 Internet Distribution of Affine Models

The networked distribution of affine models requires a two-part query—response format.
In the first part, agent clients request from lower—level server agents, models of systems
valid around specific output operating conditions. In the second part, server agents
provide the models and the input operating conditions required to operate the system

at the desired outputs.

 

Client
Query of a system System model and
model valid around required input
5’0 ii

0

 

 

 

 

 

 

 

 

 

Tier 1
Query of
Query 0f model 2
model 1 valid around
valid around
Model of Yo,2 Model
y 0,1 component 1, Component 2,
and required and required
in ut u ' ut u
‘ 7 p 0,1 ‘ ' mp 0,2

 

 

   

Figure 3.4. Query—Response Model Distribution Format

36

In a two—level network example (Figure 3.4), the User client makes a request to
the Tier 1’s server agent for a model of a system valid around the desired output
Operating conditions 370. The assembly agent determines that two model components,
the first Operating around the output y0,1 and the second Operating around 3102,
are required. The Tier 1’s client requests a model to each Of the two component
agents. Two responses are generated. In the first, the server 1 returns the model
component 1 and the required input “0,1 to Operate the first component around 90,1-
In the second, server 2 returns the model component 2 and the required input "0,2
to operate this second component around 310,2 The assembly agent uses component
model information to execute the assembly. The assembly agent also returns to the
user client, the system model and the required {'10 inputs to operate the system around

the desired outputs yo.

3.3 The Standard Model Formats

This section defines the standard dynamic model formats used in the MMM process
and the standard equations used to transform deviation variables into physical vari-
ables. In the proposed networked environment, component models are distributed as
linear models defined in deviation variables but assembled as affine models defined in
physical variables. Initially the linear format is presented, then the required variable
transform equations are defined.

The standard MMM model format can be Obtained by substituting the change Of

variables
‘ — Qi— — Jdi V < '<
y—y(t) —y0, , —— ,, (1 _ z _ n) (3.20)
dt dt
_ djfi dju -
u=ut—u, —=—,V1§ 5m
() 0 dt, dt ( J ) (3.21)

mm the affine model (3.19). This process yields a linear differential equation in the

deviation variables y(t) and ii(t),
Nay + - - - + N, (dny/dt") + Mon + - - - + Mm (ci’nn/dtm) = 0 (3.22)
Applying the Laplace Transform to (3.22) yields,
[ O n — O m _
Nos +---+an ]Y(s)+[MOs +~~+Mns ]U(s)=0 (3.23)

where Y(s) and U(s) are the Laplace transform of the deviation variables y(t) and

u(t). Defining the two polynomial matrices,
0 n
N(s) = [N03 + . - - + an ] (3.24a)
0 m
M(s) = — [Mos + - - - + Mns ] (3.24b)
and substituting them into (3.23), yields the general form
N(s)Y(s) = M(s)U(s) (3.25)

where N (s) and M(s) are (r x r) matrices.
Two external model representations are derived from (3.25). These are the dy-

namic stiffness [20] and the Transfer Function representation [21]. The Dynamic

representation
P(s)Y(s) = 0(3) (3.26)
uses the (r x r) matrix
we -—— [M(s)]‘1N(s) (3.27)
The Dynamic representation
Y(s) = G(s)U(s) (3.28)
uses the (r x r) matrix
C(s) = [New1 M(s) (3.29)

38

Affine model representations (3.26) and (3.28) are respectively analogous to linear
model representations (2.5) and (2.7). As it was shown in the last chapter for the linear
assembly process, model representation (3.26) is a convenient format for assembly
affine physical system models but it is not apprOpriate for simulation since is not
possible to solve directly for the outputs. Model representation (3.28) is used for
simulations of affine models because it can be used to compute the output given the
inputs.

The required change Of variables is obtained by defining the unit step u(t) and
substituting the equivalent time functions y; = you(t) and u; = u0u(t) into (3.20)
and (3.21),

W) = y(t) - 3'36) (330a)
u(t) = u(t) — u3(t) (3.30b)

The Laplace Transform Of (3.30a) and (3.30b) is,
Y(s) = Y(s) — y0(1/s) (3.31a)

U(s) = U(s) — 110(1/3) (3.31b)

where y0 and u0 are respectively the output and input operating vectors of the com-
ponent, Y(s) and U(s) are the Laplace transform of the output and input deviation
variables, and Y(s) and U(s) are the Laplace transform of the output and input phys-
ical variables. Equations (3.31a) and (3.31b) are the component standard equations

used to change deviation variables into physical variables.

3.4 Physical Assembly Constraints

Affine component models are specified with respect to the output and the input
component operating points y0 and uo. The resultant assembled affine model is

specified with respect tO the output and the input assembly operating point vectors

39

yo and {'10. Since component Operating point vectors y0 and uO are a subset Of the
component variables y c and uc and the assembly Operating point vectors yo and fro

equations (2.15) and (2.18) are

are a subset of the assembly variables y a and ua ,

also satisfied by the component and assembly Operating points

y0 = syo (3.32a)

T
S u0 = r10 (3.32b)

where y0 and u0 are (r x 1) and yo and 130 are (l x 1). Equations (3.32a) and (3.32b)
provide a method Of computing the output operating points of the components from
the output Operating points Of the assembly.

The process of determining system’s Operating point (yo, 1'10) follows a standard
network procedure (Figure 3.4). The client specifies the Operating point output yo
in a request to the assembly agent for a model. The assembly agent knows the
assembly constraints S and uses (3.32a) to compute the operating point outputs yo
for all the assembly’s components. These component outputs are sent to component
agents as part Of a request to each component for a model. At every level, each
component responds to a model request with that component’s linearized model and
that component’s Operating point inputs uO. The assembly agent determines the

assembly’s Operating points inputs {'10 using the connection constraint S and (3.32b).

3.5 Assembly Of Affine Physical Models

The MMM uses a systematic process to assemble dynamic matrix-based models.
These models have port pairs standardized through the two concepts that charac-
terized physical systems. The assembly process for affine models includes three steps.
These steps are: 1) Generation of the unconstrained component model, 2) Generation

of constraint equations, and 3) Generation of the assembly model.

40

The unconstrained component model is the first step. This process formulates
a diagonal matrix of component dynamic matrices. An assembly of A: component

models with a total number of r ports yields,

P1(s) O 0
Pc(s) = O '-. O (3.33)
0 Pk (s)

where P,(s) is the (r,- x ri) dynamic matrix Of the 2th component and Pc(s) is the (r x r)
unconstrained dynamic matrix of components. The unconstrained component model
relates the component deviation output vector Yc(s) to the deviation component

input vector [70(3) in the form,
Pc(s)Yc(s) = [70(3) (3.34)

The required equations to transform deviation variables into physical variables are
also available.

Yc(8) = Yc(8) - y0(1/ 8) (335a)

06(3) = Uc(s) — u0(1/s) (3.35b)
The formulation Of the dynamic constraint equations (2.19a) and (2.19b) is the sec-
ond step. The objective is to generate the constrain matrix S that relates component
output variables and assembly output variables. This relation is established by con-
sidering that output variables are constrained to be equal to other outputs variables
when ports are connected. The process is illustrated in the next section.

The generation Of the assembled model is the third step. Initially the equation

(3.35a) and (3.35b) are substituted into (3.34)
Pc(s) [Yc(3) — you/e] = Uc(s) -— u0(1/8) (3-36)

Multiplying both sides of (3.36) by ST yields,
sTPce) [Ye-(8) — yon/5)] = STUse) — sTuOU/s) (3.37)

41

Substituting (2.19a),(2.19b), (3.32a) and (3.32b) into (3.37) yields,
sTPC<s> isms) — WOO/8)] = use) — tron/s) (3.38)
Factoring matrix S in the right side and rewriting (3.38)
Pats) M(s) -— you/3)] = U.(s) — 1100/5) (3.39)
where the assembly dynamic stiffness matrix is
Pa(s) = STPc(s)S (3.40)
Defining a new set Of (l X 1) deviation variable vectors,

ms) = Yes) — 300/8) (3.41a)

Ua(s) = Uc(s) — {10(1 /3) (3.41b)

and finally, substituting (3.41a) and (3.41b) into (3.39) yields the assembly model in

deviation assembly variables,
Pa(s)Ya(s) = Ua(s) (3.42)

The assembled model (3.42) is in the MMM assembly standard format (3.26). The
process is recursive. Model (3.42) can be assembled to other standard higher order
models using the same algorithm. As in the linear case, model (3.42) must be inverted
tO Obtain the transfer function model in deviation variables. The transfer function

model is used in the simulation procedure

3.6 Example

The assembly of a mechanic transmission model and an electric generator model
is developed in this section (Figure 3.5). The electric generator has two ports. The

generator first port variables are the input electrical charge q(t) and the output voltage

42

potential e(t). The generator second port variables are the input rotational torque
71 (t) and the output angular displacement 61(t). The mechanical transmission has
two ports. The transmission first port variables are the input rotational torque 72 (t)
and the output angular displacement 02 (t) The transmission second port variawa
are the rotational torque r3(t) and the output angular displacement 03(t).

The assembly uses seven ports and three joins. Three of these seven ports are
assembly ports. The other four are component ports. The left joint connects the
charge—potential assembly port (q* (t), 6* (t)) to the component port pair (q(t), e(t)) .
The right joint connects the torque—angular displacement assembly port (r; (t), 193 (t))
to the component port (r3(t),63(t)) . The middle joint connects the component
ports (r1 (t), 01(t)) and (72(t), 02(0) to the torque-angular displacement assembly port
(ra(t),9a(t))- .

 

  

  

4‘0) I f em to). 1 9.0) mo | + arm
em 652 ‘12.")
L; ----- (— -')
9
M 7.0) 120) 30L»

 

 

13(1)

 

Figure 3.5. Tfansmission—Electric Generator Assembly

The external nonlinear models for the generator and the transmission have respec—

tively the general form,
f9 =(e,61,é,01,... i41,71.41.+1,---) = 0 (3.4321)
fr = (92,93,92937'“ ,T2iTsaT'2sT'3w“) = 0 (3.43b)
where fg(t) and f7(t) are vectors of functions. The user requests a model of the

assembly by defining its output operating points. These are the operating voltage

63, and the operating angular displacements 03‘“, and 9"; 3. Using this information,

43

the assembly agent computes the generator output Operating points e0 = e; and
00’, = 60,0, and the transmission output Operating points 602 = 60,0 and 60,3 2 93,3.
Component agents have at this point all the information to provide the component

models and the physical variable transformation equations. For the generator,

 

 

1311(8) P12(3) _ é1(5) _ T109) 3

i 1013(8) 1914(5) l _ i 3(5) l _ i Q(3) l (3.44 )
é1(8) 1 = , 91(3) _ 90,1 ]l

l 1%) 1 _ E(s) l l e. s (3'44”

27:1(3) 1 = ' T1(s) _ 70,1]1 c

l 62(8) . _ 62(5) l i «1., s (3'44)

where (2(3), T1(s), E(s) and 91(3) are respectively the Laplace Transform of the
deviation time variables q(t), r1 (t), e(t) and 6] (t). The complex variables Q(s), T1(s),
E (s) and 61(5) are respectively the Laplace Transform of the physical variables q(t),
71(t), e(t) and 01(t). The complex variable functions 1311(3), 1912(3), P13(s) and
P14(s) are the dynamic elements of the generator and finally the constants r0,1 and
q0 are respectivelly the Operating input torque and the Operating input charge.

The mechanical transmission model is,

 

P210?) 1322(8) : E:92(S) : 7:2(8) a

i 1323(8) 1024(3) l i 63(3) l i T3(5) i (3.45 )
62(3) ‘ __ _ e2(8) _ 90,2 _1_ 5

i (93(3) _ _ . 93(8) l [ 90,3 3 (3.4 b)
73(3) - _ ' T2(3) __ 70,2 l c

[T3693 . — _T3(8) l i T0,3 i 5 (3.45 )

 

 

where 92(3), 73(3), (93(3) and T3(s) are respectively the Laplace Transform of the
deviation time variables 32(t), r20), 63(t) and 23(5). The complex variables Q(s),
T1(s), E (s) and 01(3) are respectively the Laplace Transform Of the physical variables
q(t), “q(t), e(t) and 91(t). The complex variable functions P21(s), 1922(3), P23(S) and
P24(s) are the dynamic elements of the transmission and finally the constants r0,2

and r0,3 are the Operating torque inputs.

44

The first step is to generate the unconstrained component model. Initially the

unconstrained component matrix is generated. This matrix is,

1911(3) 1312(8) 0 0
_ 1013(5) 1314(8) 0 0 .
136(8) _ 0 0 1021(8) 1022(8) (346)
0 0 1323(8) 1324(3)
The unconstrained component model is,
P11(8) 1312(8) 0 0 é_1(3) C13(8)
P13(3) P14(3) 0 0 _E(S) = 9(8) (3 47)
0 0 1321(8) 1322(8) €209) 732(8) '
0 0 1323(8) 1324(3) 93(8) T3(8)

The second step is to generate the constraint assembly equations. Initially the
time dependent equations that relate component output to assembly outputs are

written. These equations are,

6t=() 6%)
91(t)=92(t9)= a(t) (3-48)
93t=() 9305)

Additionally, conservation of energy at each joint requires

q(t)€(t) = (1*(t)€*(t)
Ta(t)9a(t) = T1 (109105) + 720092“) (3-49)
73(t)93(t) = Té‘ (095“)
From (3.48) and (3.49), follows that
(1“) = C1%)
u(t) = 7'1“) + T205) (3-50)
730‘) = 7;? (t)

Applying the Laplace Transform to (3.48) and (3.50) results,

e1(5) E*(8)
E(s) = s 90(5) (3.51a)
92(5) 6*(8)
93(8) 3
T1(5) * 8
ST Q“) 3(3)) (3.51b)
T2(3) Tics)
73(8) 3

45

where the constraint matrix S =

and 7*

a:
where qO, To a 0 3

 

Cot-‘0
CHOP“

 

T0,2

 

7"0,1

HOOD

 

- 70,3 _

. Similarly, for the operating points,

*

7-O,a

T0,3

are the operating inputs of the assembly.

(3.52a)

(3.52b)

In the third step, the assembly model is generated. Substituting (3.44b), (3.44c),

(3.45b) and (3.45c) into (3.47) yields,

 

 

:{ 91(3) - 00,1 I: 3‘
' E(s) e 1 - '
P 0 — =
ml. 92(3) 00,, 3 l i
I 93(3) 603 _ J I
Multiplying both sides of (3.53) by ST yields,
:{ 91(3) - 60,1 1 I. I 73(5)
' E(s) e I ' ' T Q(3)
STP 0 _ =
Ml. 92(3) 90,, 3 l is T2(3)
l 93(8) _ 90,3 - J. I 73(8)

 

 

T1(5)
Q(S)
T2(S)
T3(5)

 

_ 70,3 J

 

70,1
90
T0,2

 

Replacing (3.51a), (3.51b), (3.52a) and (3.52b) into (3.54) yields,

T E*(s) 63 ll .‘ Q*(s)
SPC(3 )(s 90(3) -s 60,, ;}=( Ta(3)
l 93( s) 933 J l T518)
Factorlzmg the matrlx S 1n the r1ght Side of (3 53) yields
E*<s> z; 1). .r em
Pa(s){ ea<s) o... ;}=( Ta(s)
l 93(8) 033 J l 73(5)

qo

TO,a
*

70,3

T0,1

a: | ._.

 

0: I r—t

a:

Co I t—I

CI) I t—l

(3.53)

(3.54)

(3.55)

(3.56)

where the assembly matrix Pa(3) = STPC(3)S. Defining a new set of variables ,

 

{73*(3) E*(3) - e; 1
9a(3) = 9a(3) — 90,0 ; (3.57a)
93(3) 93(3) . 93,3
62*(8) ms) ' q; ,
_a(3) = Ta(3) — 70,0 ; (3.57b)
T§<s> ms) . 73:3
and replacing them into equation (3.56) yields,
Plus) 1212(5) o E(s) 62*(s)
1913(8) P14(S) + 1721(8) 1022(8) (_9a(8) = 22(8) (3-58)
0 1023(5) 1024(8) 93(8) T; (8)

Equation (3.58) is the the assembly of the generator model and the mechanical trans-
mission model. This resultant model is in the standard format defined in (3.26).
The variable transformation equations (3.57a) and (3.57b) are in the standard format
(3.31a) and (3.31b). The process is recursive and the model (3.58) can be assembled

to other standard higher order models using the same algorithm.

3.7 Summary

In this chapter, two important contributions to mechanical engineering are presented.
The first contribution is a method to assemble physical port-based affine ODES around
an equilibrium operating point. Affine systems often result of local linearization about
an operating point. In this case, the local system model is linear in deviation vari-
ables and non-linear in physical variables. Because assembly constraints are always
given in terms of physical variables, an explicit method to apply physical assembly
constraints for affine physical models was deveIOped. This method addresses one of
the most common nonlinear systems in mechanical engineering. Using this method
many practical physical nonlinear system models can be obtained by assembling the

affine approximations of their subsystems models.

47

The second contribution is an explicit method to obtain subsystems operating
points from the system assembly operating point outputs. This method for the first
time provides an explicit, closed form solution to the general operating point problem.
In a recursive closed form, the method specifies information required to roll-down an
operating point specification through every subsystem to all the lowest level compo-
nents of a system. For each lowest level component, the method provides the exact
port output values where the operating point representation of the nonlinear com-
ponent model should be developed. The method then provides an explicit method
for computing operating point inputs of the component level and use these inputs to

compute the operating point inputs of the assembled system model.

48

CHAPTER 4

Volterra Model Formats

4. 1 Background

Physical nonlinear dynamic systems have traditionally been modeled using nonlin-
ear ordinary differential equations (ODEs). One solution technique for single-input,
single-output (SISO) nonlinear ODES is obtained through a Volterra transfer func-
tion. Applying this solution technique for multi-input, multi—output (MIMO) non-
linear ODEs is a complicated issue due to the difficulty of analytically determining
MIMO Volterra models.

MIMO Volterra transfer function models are multi input—output representations
used to describe nonlinear behavior. A Volterra model is nonlinear with respect to its
input but linear with respect to its parameters [24]. Volterra models have been used
to obtain solutions for SISO nonlinear systems in areas such as dynamic behavior of
offshore structures [25], power amplifier behavior [14], rheology [15], nonlinear model
reduction [16] and others. Areas such as sub-harmonic nonlinear behavior [26] and
noise characterization [27] use Volterra models for describing the behavior of multi-

input single-output nonlinear dynamic systems.

Although an analytical procedure to obtain Volterra transfer functions for SISO
nonlinear ODEs is available [28], [10], no general analytical procedure to obtain a

MIMO Volterra transfer functions from port—based, external nonlinear ODEs has

49

been published. This is a complicated issue when the number of system’s inputs and
system’s outputs are not equal. The orientation of existing works in MIMO Volterra
models is to experimentally identify the model kernels [29], to approximate weakly
nonlinear MIMO systems [30] and to determine error bounds [31].

In this chapter, a procedure to analytically obtain MIMO Volterra models from
port-based nonlinear ODEs is derived. This procedure provides a solution to this
class of MIMO nonlinear ODEs. First, the standard nonlinear port-based ODEs used
in this procedure are described. Second, the Volterra MIMO models are prasented.
Third, a nonlinear operator required to obtain the MIMO Volterra models is defined.
Fourth, the procedure to obtain MIMO Volterra model expansions from port-based
nonlinear ODEs is presented. Fifth, an example of a two-port nonlinear ODE is shown.
The Volterra model response is compared to the nonlinear ODE model response using

a Simulink simulation.

4.2 Port Based Nonlinear ODE Models

A port—based nonlinear ODE model has equal number of inputs and outputs because
an input—output pair constitutes each port. This section considers nonlinear port-
based models in the form,

. dm - dm
f1(y17y17"'Wail7"'iyrvyra"°—az7%t) =U1
- (4.1)

fr (yty'hm%,---,yr,yr,...‘13-?) = Ur

where the number of ports 1 S i g 7‘; f,(-) is the 72th port’s nonlinear operator; u,

is the ith input and y,- is the ith output of the system. , Many physical models of

nonlinear engineering systems take this form when modeled using the energy-based
Lagrange approach [32].

An example is a two—port system (Figure 4.1) composed of two masses, one nonlin-

ear spring, one linear spring and one linear damper. The mass m1 is connected to a

50

fixed point through the nonlinear spring with force f N L = 161(3)1 + yi). The mass m2
is connected to the mass ml with a linear damper having damping coefficient b1 and
a linear spring with constant k2. The first port uses the input-output pair (u1,y1)
where u1 is the input force applied to the mass and y1 is its output displacement. The
second port uses the input-output pair (212,312) where 112 is the input force applied to

the mass m2 and y2 is its output displacement.

 

 

“10) 712(1)

—> —>

k2

Nonlinear _ W m
Spring - '_ 2
b. I_,

MO) 150)

 

 

 

 

 

 

 

 

 

 

Figure 4.1. Double-Mass Spring Damper-Nonlinear System

The nonlinear, port—based, ODE model for the system (Figure 5.3) in the required
format (4.1) is

"“1171 + b1(91— 92) + k1y1+ (“11/13 + k2(y1 ‘ yz) = "1

.. . . 4.2
m2y2+b1(y2—yl)+k2(y2—y1)=u2 ( )

4.3 Volterra Models

In this section two nonlinear model formats are presented. The first model format
is the MIMO Volterra transfer function representation. This model representation is
used in the simulation of port based nonlinear dynamic models. The second model
format is the MIMO Volterra dynamic model. This model representation is used in
the nonlinear MMM assembly of port based dynamic physical systems.

A Volterra transfer function model is an extension of the linear convolution integral

approach for time—invariant systems. Any nonlinear time-invariant system can be

51

represented as an infinite sum of multidimensional convolution integrals of increasing

order [10]. For a SISO system this is
00
y(t) = 491(Tl)u(t — T1)d7'1+

[0 [092 (71,72)u(t — 7'1)u(t — T2)d7’1d7'2+ (4.3)

000000

f0/0f093(71,Tgr3)u(t—rl)u(t—72)u(t—T3)d71d72d73+,,,

where y(t) represents the nonlinear system output, u(t — T1) is the input associated
to the ith dimensional convolution integral, and 9,-(7'1, . . . , Ti) is the ith kernel of the
system. All the higher order kernels can be analytically derived from the first order
linear kernel 91(71) [28]. This kernel is the traditional linear unit impulse response
of the system. The second order kernel, 92(71, 7'2), is a two-dimensional function of
time and is the response of the system to two independent unit impulses applied at
two different points in time. The kernel 92(71, 72) is a function of both time and one
time lag (71 — 72) . The kernel 93(71, 72, T3) is a three-dimensional function of time
and is the response of the system to three independent unit impulses applied at three
different points in time. The kernel 93(71, 72, T3) is a function of both time and the
two time lags (71 — 72) and (T2 - 7'3) [33].

The ith integral in (4.3) is called a ”degree-2' homogeneous system” and its respec-

tive output

0000

312' =// 91(71w” ,Ti)u(t—71l°“u(t-Tild7‘1"'de' (4-4)
0

0
is called the ”degree-i homogeneous output” (4.4). Notice that if a new input au
where a is a scalar, is applied, the resultant output is g], = anyi [10]. The SISO
response y(t) in (4.4) is then a sum of partial outputs provided by different ”degree-2'

homo eneous s stems” in the form
')

y(t) = y1(t1) + 31201. t2) + 31361, t2, t3) + ' " (4-5)

52

where,

00
91 =/0 91(71)U(t—71)d71
00 oo

92 *4/0/ 92(71,T2)U(t—T1)u(t-T2)dTldT2

oo oo 03
332/0 /0 [0g3(Tl,T2,T3)u(t—’rl)u(t—’r2)u(t—'r3)d7'1d'r2dr3

and gi('rl, - - - Ti) is the ith' kernel.

The Laplace transform of a multi-time variable system is an extension of the single-
time Laplace transform. Given an ith dimensional function of time 9,-(7'1, - - - Ti) that

is one—sided in each variable (0 S T,- S 00) , its Laplace transform is defined as [10]
00 oo
Gi(sla°” 18i)=/ 0”] 92.017.” ,E)e_slT17'” ie—siflngla°” 1dT1’ (4'6)
0 0

where (1 _<_ k S 2') the set of complex variables 3k = 0k + wkj. Applying the Laplace

transform to each term in (4.5) yields,

Y1(31) = 01(81)U(31)
Y2(31,32) = G2(81, 82)U(81)U(82)

Y3<319 82, 83) = 03(31, 82, S3)U(31)U(32)U(33) (4.7)

where U (3,) is the Laplace transform of the input u(ti) and Gi(31,- -- ,31.) is the
multi-complex variable Laplace transform of the ith order kernel. The SISO Volterra
transfer function is obtained by summing all the rows in (4.7)

Y(81,82, ' ' ') = 01(81)U(81) + 02(81, 82)U(81)U(52)

(4.8)
+G3(31, 82, S3)U(31)U(32)U(53) + . ..

The model representation (4.8) is used to simulate the response of SISO nonlinear

systems [10].

53

The MIMO extension of (4.7) for a system with 7" ports is,
Y1(81) = G1 (U(81ll

Y2(81, S2) = G2 (U(31)U(82))

(4.9)
Y3(81, .92, 83) = G3 (U(31lU(82)U(53))

Each element Gi(o) in (4.9) is a functional that operates on a set of (7‘ X 1) input
vectors, yielding V(£ i S 00) the (r x 1) ”degree—2' homogeneous output” vector,

T
Yi(517”' 152:): [K,1(813”' 137:)1'” 13/137481)”. 181)] (410)

The Volterra transfer function for a port-based system with 7" ports is an extension

of (4.8) and is obtained by adding all the (23(0) functionals in (4.9). The response is

a vector in the form,

Y(81,-°')= G1(U(31))+ G2 (U(81), U(S2)) + G3 (U(sl), U(32), U(33)) + . --
(4.11)
where U(31),U(32), - -- , are (r x 1) input vectors and Y(31,-~) is the resultant

(r x 1) output vector.

The MIMO notation in (4.11) can be simplified if the following input argument

list is defined,
U1 '3 U(Sll

U2 3 {U(51),U(82)}

A (4.12)
U3 : {U(81)1U(52)1U(33)}

Using (4.12) into (4.11) yields,
Y(81,82, - ' °) = G1(U1) + G2(U2) + G3(U3) + ° " (4.13)

In the same form, the output notation is further simplified using
Y(1) 3‘ Y1(31) = G1(U1)
Y(2) é Y2(81, 82) = G2(U2)

A (4.14)
Y(3) = Y3(51, 52, 83) = G3(U3)

54

Using (4.14), the MIMO Volterra transfer function model in (4.13) can be written as,
Y(31,--') =Y(1) +Y(2) +Y(3) +°°- (4.15)

This model format is used in the simulation of nonlinear port based physical systems.

A procedure generates each of the Y(i) [34]. Each Y(,-) has an unique standard
format that can be represented in function of the first order linear Operator G1(U1) =
G1,(1)U1. The term G1,“) is the (r x 7‘) transfer function matrix Obtained from the
linearization of the system nonlinear ODE in (4.1) around an equilibrium Operating

point Op. The matrix G1,“) can also be Obtained by inverting the dynamic matrix

m=1<:—::)11<—>1[Hit—>1

where (1 S i, j S 7‘) and each element in the brackets is a (r x r) matrix Of constants.
The sub-index (1) in G1,“) or P1,“) indicates that these matrices are function of
a single complex variable 31 For any sufficiently differentiable ODE, the first three

standard Volterra Operators are [34],

Y(l) = G1,(1)U1 (4.178.)

Y(2) = ‘G1.<2)P2,(2)¢’ (GMDULGMIWI) (41713)

P3,(3)¢(G1,(1)U1’G1,(1)U1’G1,(1)U1)+
P2,(3)¢lG1,(1)Ula G1,(2)P2.(2)'(’(G1.(1>U11G1,(1)U1)l

where GM”) is the result of inverting the dynamic matrix P110), and 13(0) is a non-

(4.17c)

 

Y(B) = 43143)]

linear Operator. Any dynamic matrix Pu,(v), includes 2) complex variables, 31, - - - ,sv

and can be derived from (4.1) using,

I 6‘1“ 3.5.. q
: [WW [(3)1343

Pu:(’U) : a 4 617, 2 ] a V, 3
'1 [(554)].Opsm + _(W)]10,.(v)+...

J .7

(4.18)

 

where 3(1)) = 31 + - - - + 31,.

55

The Volterra Dynamic Model is other Volterra format. The Volterra dynamic
model is composed by a set Of multi-complex variable dynamic matrices. For the ith

component model that includes q expansion terms, this model has the form,

{Pi,1,(v)1 P232411), ' " ’Pi,q,(v)} ' (4-19)

Model (4.19) is used for the MMM to assemble nonlinear physical models. For any
Pi,u,(v) in (4.19), the index 2', represents the component, the index u represent the
order Of the expansion and the index 12 represents the number of complex independent

variables used. Each Pi,u,(v) matrix satisfies,

Pi,u+1,(v) = me. (Pam) (4.20)

where each element of the (r x 7‘) matrix Dy,y,y.-" (Pi,u,(v)) is the first partial deriva-
tive respect the vectors y, y, y, - - - Of its corresponding elements in the (r x r) matrix
Pi,u,(v)' Equation (4.20) is derived from (4.18)by rewritten it in general form V(1 S u).

In example, for the ith component, the dynamic matrices,

Pi,2,(v) = Dy,y,yr" (P231410)
13,33)”, = DEW... (P,,2,(,,))
P1,4,(11) = Dy,y,yr" (13233160)

The Volterra series is an infinite power series with memory and suffers from the
problem Of a limited zone Of convergence. For the Volterra representation to prop-
erly describe a subsystem system model, the Volterra series must converge over the
input/output variable range used in the assembled system. This issue is recognized
in multiple works: Frank [39], Brilliant [11], Christensen [40], Boyd [13], Czarniak
[41], Tomilson [43], Chatterjee, [42] and others. The series convergence depends on
system parameters [42], BIBO stability [40], amplitude Of the input applied [13], trun—
cation point of the series and other factors. One necessary condition for existence of

the series and its convergence is that the linear term cannot be zero [39]. Boyd [13]

56

showed that the radius Of convergence of a Volterra series is directly related to the
amplitude of the input applied. This radius Of convergence can be presented in the
form |u(t)| < p, where u(t) is the input vector applied tO the system and p is the
radius Of convergence. Sandberg [12] has also shown that a truncated Volterra series
provides a uniform approximation to the in finite Volterra series on a ball of bounded
inputs for a large class Of systems. SISO Volterra series have been recognized as
a useful approach tO engineering modeling. This work assumes a useful convergent

Volterra model exist and converges over the domain Of model application.

4.4 The Nonlinear Multiplicative Operator

A nonlinear Operator 10(0) is used to build the Volterra functionals in (4.13). This
Operator uses as argument an arbitrary set of vectors (r X 1) in the form,

U1 U1 w1 21
R: u: ; ,v= 3 w: ; ,---,z= ; (4.21)

24,— vr wr Zr

to perform the Operation rb(R) that yields the (r x 1) vector

(U1) ' (’01) ' (101) ' ° ° ' (31)
13 = ME) = 2 (4.22)

(Ur) ' (Ur) ' (107‘) ° ’ ' ' (27')
Each element in the vector p is a scalar formed by the product Of corresponding
elements Of the argument vectors u,v,w, - -- ,z. Restated, the ith element Of the
vector 13 results of multiplying the ithelements of the vectors u, v,w, - - - ,z , that is
151: (Ui)-(v4)-(w4)---(v4)-

The result Of Operating on a set of (r x 1) degree-i homogeneous output vectors by
the multiplicative Operator 10(0) is defined to have the standard notation TWA", ,w).
This process also satisfies (4.22), but the resultant vector notation includes the ad—
ditional sub—index (a, b, - -- ,w). This sub-index is used to identify the number of

argument vectors with equal degree—2' homogeneous outputs. In the notation, a is

57

the number Of degree—1 homogeneous output argument vectors, b is the number of
degree-2 homogeneous output argument vectors and so on. Finally w is the number

Of argument vectors having the maximum homogeneous degree.

Two examples are presented to illustrate the multiplicative Operator ¢(o) notation.
In the first example We) Operates on the combination Of two (r x 1) degree-1 homoge—
neous output vectors. In the second example 13(0) operates on the combination of one
(r x 1) degree-1 homogeneous output vector and one (1' x 1) degree-2 homogeneous
output vector. These examples will illustrate the use of the required sub-index when

different combination of degree are present in the operator’s argument list.

In the first example, the (r x 1) vector

~

Y(2) = 14')(Y(1),Y(1)) = [(3/(1),1) “(Y(1),1),“' 1(Y(1),r) '(Y(1),r)lT (4-23)

is the result Of Operating two degree-1 homogeneous (r x 1) output vectors Y(l) =

“Y(UJ)’ - - - ,(Y(1),T)]T] (Figure 4.2). Each element of the vector \7 is Obtained by

(2)
performing \7’(1 S 2' S r) the Operation (YUM) - (Y(llfl')‘ The sub-index (2) in ?(2)
indicates that it is the result Of Operating on (2) degree—1 homogeneous output vectors.
Because the indices in parentheses end after the first entry, no higher degree terms
are indicated.

The resultant output Y can be also presented as a function Of the external

(2)
input U1. Substituting the first row Of equation (4.14) into (4.23) yields,

Ya) = 4 (G1<U1).G1<U1)> (4.24)

Using (4.24), Volterra functionals (4.14) that includes Y(Q), can be written as a

function of the system external inputs.

In the second example, the (r x 1) vector

Y(1,1) = ¢(Y(1),Y(2)) = [(3/(1),1) ' (Y(2),1), - " 1(Y(1),r) ' (Y(z),r)lT (4-25)

58

 

 

 

 

 

 

 

 

 

 

 

i; E Y?” ‘ A
a wa ' : (11(1),!);(Y(I)J) Y(2)
' E
1: Ya)» l

 

 

 

 

 

 

 

 

 

Figure 4.2. Block Diagram Representing the Operator 1,0(Y(1),Y(1))

is the result of Operating on a. combination of one degree-1 homogeneous (r x 1)
output vector, Y0) = [(YmJ), - - - , (Y(l),,)]T with one degree-2 homogeneous (r x 1)
output vector, Y(Z) = [(I/(z),1),--- ,(1/(2),,.)]T] (Figure 4.3). Each element of the
vector 170,1) is obtained by performing V(1 S i S r) the Operation (YUM) - (Y(Z),i)-
The sub—index (1,1) in Y(LI) indicates that it is the result of Operating on one degree-
1 homogeneous vector and one degree—2 homogeneous vector. Because no additional

indices are included beyond the first two, no higher degree arguments are present.

 

 

 

 

 

 

 

 

 

 

 

Y(l) é Y(l),1
E Y(l),r : (Ylll,l):(Y(2).1) film)
i E
#32 Y” I (3)3132»)
c: Y(2).r

 

 

 

 

 

 

 

Figure 4.3. Block Diagram Representing the Operator 111(Y(1),Y(2))

The resultant output Y(1,1) can be also presented as a function Of the external

59

input U1 and the input set U2 by substituting the first two rows Of equation (4.14)
into (4.25) to yield,

Y(1,1) = 1!) (G1(U1), G2(U2)) (426)
Using (4.24) and (4.26), Volterra models including Y(lJ) can be written as a function
of the system external inputs. This capability is important because the Volterra

models Obtained from MIMO port-based ODEs will include (4.23), (4.24) and other

similar higher order terms.

4.5 Generating Volterra Models From Port-Based ODEs

The Volterra model Of a port-based MIMO nonlinear ODE model in the form (4.1),
is Obtained in this section. Assuming that this system is Operating around an equilib-

rium point (op) located at the origin (yo, no) = (0,0) and using the vector notation
a" = [on ,an]T (4.27)

a Taylor series expansion is applied around this equilibrium point V(1 S 2', j S r) and

V(1 S k S m), (m is the maximum time derivative of the outputs) to yield,

iWSEOIOM[(3%)le[(3)14lv”+-~>l=u

”:1 J J .7

Defining the (r x 1‘) matrix of constants,

a";
[big] (wk) 2 KW lop)] (4.29)

For an specific k and 77., each constant bi j in the matrix, is the nth partial derivative
Of the scalar function fi (y, y, - - ) respect the variable [(1ky;I /dtk] evaluated at (op).

Using a constant a, a new input in the form

fi=au=a[ul(t1),--- ,u,.(t1)]T (4.30)

60

is defined. If this new input is applied to (4.28), the resultant output 5: is a MIMO

extention of the SISO form shown in [28],

5’ = @101) + 023’2(t1,t2) + O3Y3(t1,t2,t3) + ° ' ' (4-31)

Equation (4.31) is a direct result Of the ”degree—2' homogeneous system” properties
(4.5) and is extended here to the vector case. Defining y(i) = yi(t1, - -- ,ti), (4.31)
can be written in the form,

5' = ayu) + (12y(2) + a3y(3) + - - - (4.32)
Substituting (4.30) and (4.32) V(l S z',j S 7‘) into (4.28) yields,

lbi,jl(1’0)(a)ll)+azy(2) + ' ' ) + lbi,j](1,1)(a5fl) + (125(2) + ' ' > + ° ' ° +

2 2
21! lbi,jl(2,0)(ay(1)+a2y(2)+ ' ' ) + 211. lbi,jl(2,1) (”if 025(2)+ ' ° ) + ' "+
3 . . 3
31—] [bi‘j](3’0)(ay(1) + a2y(2) + - - -) +31] [bi,j](3,1)(a3[1) + 03y(2) + - - ) + = on

(4.33)

Terms multiplied by a particular ap, V(1 S p S 00) in (4.33), are independent from all
others because a is an arbitrary constant. This property is used by Schetzen [28] to
determine all the SISO Volterra transfer functions in (4.8). Extending this procedure
to the MIMO case, all the Volterra Operator Gp(o) in (4.11) V(1 S p S 00) can be
derived from the terms that are multiplied by the constant up.

The first term in the expansion (4.11) is the linear Volterra Operator G1(o). This

term is derived from the elements multiplied by 01,

m
' " d y“) _ 4 34
[bi’leOPfUJr[bi’il(1.1)’ll)+[bi’jlmfilnl’ "'+lbz',jl(1,m)"&2"r' ‘ ‘1 ( ' )
where m is the maximum output derivative. Define the time dependent I“ order

vector Operator,

p1(x) = [bi,jl(1,0)x + [bi’j](1’1)x + [bi,j](1,2)i + - - - + lbi,jl(1,m) g3? (4.35)

61

where x is an arbitrary (r x 1) vector, to write (4.34) in the form,
P1 (y(1)) = u (4.36)
Applying the Laplace transform to (4.36) with respect to t1 variable yields,
P1,(1)Y(1) 2 U1 (4.37)
where Y“) and U1 are respectively the Laplace transform of y ( 1) and u and,
131,(1) : lbi,jl(1,m)8in + ' ° ' + lbi,jl(1,2)sf + lbi,jl(1,1)31 + lbi,jl(1,o)

is a (r x r) dynamic matrix Of polynomials in the complex variable 31. The elements
between parenthesis in the sub-index of the matrix P1 (1) indicates that this matrix

is function Of one independent complex variable. All the above elements of P1 (1) in

brackets are (r x 7') matrices Of constants. Inverting P1 (1) yields,
Y0) = GMDU1 (4.38)
where
—1
2
G = b.. sm+---+ b.. s +b.. 5 +b.. (4.39)
1,0) l W](I,m) 1 l 1’3](1,2) 1 l z’J](1,1) 1 l 23]“,0)

is the transfer function matrix. Equation (4.38) is equivalent tO the linearized model
of the nonlinear system (4.28).

The second term in the expansion (4.11) is the nonlinear Operator G2(U2). This
term is found by equating the elements that are multiplied by the coefficient a2 in

both sides Of (4.33),

. dm 2
lb‘ij](1,0)bl2)+ lbiJlmfh) + ' ' ' + lbiJlUm) mm +

m 2 (4.40)
2 1 - 2 1 d yo) _

211. [bi,jl(2,0)y(1) +21 lbi,jl(2,1)y(1) + ' ' ' +21 lbz‘JlQmfifim‘ "

Using the vector notation (4.27) and the operator (4.22), the (7‘ x 1) vector
i = y 2 (4 41)
(2) (1) '

62

The definition (4.35) of Operator p1(o) and (4.41) are applied to (4.40) to yields,

p1 (y(2)) + p2 (3(2)) = 0 (4.42)

where the second order functional

m 2
p2 (x) = 71; [l2,.,].](2’0)1<2 + 21, [(21.,J.](2,1):'c2 + - - 4 + 51, [b,,j](2,m) 437%— (4.43)

Applying the Laplace transform tO (4.42) with respect to two independent time vari-

ables and using the Theorem 2.3 in Rugh [10] yields,
P1.(2)Y(2) + P2,(2)Y(2) = 0 (444)

where the (r x 7‘) matrices

PM?) = [b,,j](1,m) (31+32)m + - - - + lbw-l (31+32) + [bz.,j](1,0)

(1,1)

P242) : 21! lbml (2,m)(31+52)m + ' ' ' + 2%’ [bin] (31 + 32l+ 211' lbm'l

(2,1) (2,0)

The vector if”) in (4.45) is defined in (4.24). Substituting this result into (4.45)

yields
P1,(2)Y(2) + P2,(2)’l/) (Y(l)’Y(l)) = O (4.45)
Solving for the (r x 1) vector Y(Q) by inverting the matrix P1 (2),
Y(2) -_— —G1,(2)P2,(2)2/2 (Y(l),Y(1)) (4.46)
—1
where G139) = lP1,(2)l . Finally substituting (4.38) into (4.46) yields,
Y(Q) = ’G1,(2)P2,(2)¢(G1,(1)U11G1,(1)U1) = G2(U2) (4.47)
The output vector Y(2) is Obtained (Figure 4.4) by multiplying each of the two
inputs U1 vectors by the G1 (1) matrices. The two Y(l) resultant outputs, are Oper-

ated on by 2,0(0) to yield the vector ?(2)' This vector is then Operated by the matrix

P2’(2). The resultant output is finally Operated by the matrix —G1’(2).

63

 

 

 

 

 

 

 

Figure 4.4. Nonlinear Operator G2 (U2)

The third term in the expansion (4.13) is the nonlinear Operator G3(U3). This

term is found by equating the coefficients of a3 in both sides of (4.33),

m
d
lbi.jl(1,0)3(3)+ ' ' ' + [big] #+

(1,111) dt
2
1 1901134211 l 1 4’" 1901134211 1
+2! lbi,jl(2,0)i 3 f+' ' ' +71 lbid'l(2,m)E"r ' f
1 3101.310), J 19(1).ry(2).r J
m 3
1 3 1 d y 1 _
3'! lbi,jl(3,o)y(1)+' ' ' + 3‘! lbiij](3,m)?"$_l _ 0
The vectors
1311.134211]:
y(1,l) = | I ?
{y(i),ry(2),rJ

~ _ 3
y(3) _ yo)

and the functional (4.35) are substituted into (4.48) and,

p1 (y(3)) + W (904)) + p3 (37(3)) = 0
where the functional

7"
p3(X) = 31 lbm'lwyo)x3 + 31! lbi1jl(3y1) 5‘ + 3%! [big] (37",) (171733

64

(4.48)

(4.49)

(4.50)

(4.51)

(4.52)

Applying the multivariable Laplace transform to (4.51) with respect three indepen-

dent time variables yields,

P1,(3)Y(3) + P2’(3)Y(1,1) + P3,(3)Y(3) = 0 (4.53)

where the (r x 7‘) matrices,

P1,(3):lbml(1,m)(5<3>)m+ ' ' ' +lbidl<1,1)(3<3)) + lbiJlum
P2,(3):7lllbivj l(2,17«.)(‘("(3))m+ ' ' ' +§1llbi,jl(2,1)(s(3)) + 2-lbiJl(2,0)

P3,(3)=§1!lbi,jl(3,m)(3(3))er ' ' ' +§Ylbi,jl(3,1)(5(3)) + Elllbmlm)

_lr—I

and 5(3) = 31 + 32 + 53. The nonlinear operator (4.22) is used to evaluate the vector,
Y3 = ¢ (G1(U1), G1(U1), G1(U1)) (4-54)

Substituting (4.54) and (4.32) into (4.53) yields,
P Y(3) + P2,(3)¢ (G1(U1), G2(U2)) +

1,(3)
(4.55)
P3,(3)¢ (G1(U1), G1(U1), G1(U1)) = 0
Solving for the vector Y(3) by inverting the matrix Pl (3) yields,
Y(3)=-G1,(3)P3,(3)¢ (G1(U1), G1(U1), G1(U1)) —
(4.56)

G1,(3)P2,(3)¢ (G1(U1), G2(U2))

1
where G P . Finally by substituting (4.47) and (4.38) into (4.56) yields,

1,<3):l 1,(3>l

Y(3)= G3(U3) = “G1,(3)P3,(3)¢ (G1,(1)U1’ G1,(1)U1’ G'1,(1)U1) _

Gl,(3)P2.(3)‘/’(G1,(1)U1vG1,(2)P2,(2)¢(Gl,(1)UvGl,(1)U1)

The addition of two third order components (Figure 4.5) is required to compute the

(4.57)

output vector Y(3)' Using a similar procedure, the ith Volterra functional Gi(Uz-)
can also be calculated. Any Pu“) matrix of this Volterra operator can be calculated

using the equation,

PW) = 5.11 [bi’j](u,m) (my! + . . . + 51! [bi’jla 1) (3m) + 51! [52.0.] u, (4.58)

 

Figure 4.5. Nonlinear Operator G3 (U3)

where 3(1)) = (51 + - - - + s").

The Volterra functional in (4.13) requires an infinite number of expansions. A
realizable Volterra functional uses only a finite number of expansion terms and is
a truncated version of the general Volterra model (4.13). For example, a Volterra
functional truncated at the third expansion is obtained by adding the Volterra models
(4.38), (4.47) and (4.57)

Y1—3 = G1(U1) + G2(U2) + G3(U3) (4-59)

To obtain a truncated time solution of (4.59), the inverse Laplace transform is applied
sequentially with respect to the time variables tl,t2 and then t3. The result is the
multi-time variable truncated solution y1_3(t1, t2,t3). Finally the condition t = t1 =

t2 = t3 is applied to obtain the single variable solution y1_3 (t) [28]-[10].

66

4.6 Example

The Volterra transfer function for the nonlinear ODE model of a mass—spring system
(Figure 4.1) is obtained in this section. The two masses m1 and m2 , are given different
values to introduce asymmetry in the responses yl and 312. Because the nonlinear
spring force dominates the linear spring force, ||k1(y1+y:1”)ll 2||1c2(y2 - y1)||, the spring
parameters are selected as [C] = 8 and k2 = 4to make the nonlinear spring dominate
the system response. The damping parameter is selected as b1 = 2 to increase energy
dissipation rate and decrease settling time. Substituting these parameters in (4.2)

yields,
2371 + 23), — 292 + 8y1 + 83/13 + 4y1 - 4y2 = “1
Q2 + 292 — 2371+ 4y2 — 4y1= 142

Following the procedure (4.30) through (4.33) yields,

(4.60)

2 o .. 2.. 2 -2 . 2.
l0 1l(ay(1)+a’l2)+ )+[—2 2l(0‘y(1)+a’l2)+ )+

12 —4 (0y +02y +---) + 8 0 (ay +042)! +~--)3 =au
—4 4 (1) (2) O O (1) (2)
(4.61)
The linearized transfer function matrix G1 (1) is found by equating the elements
multiplied by a1 in (4.61),
2 O .. 2 —2 . 12 -4
l0 1l’ll)+l—2 2 l’ll) + l—4 4ly(1)" “ (4'62)
Defining the operator
2 0 .. 2 —2 . 12 —4
p1(x)—[O 1]x+[_2 2]x+[_4 4]x (4.63)
where (xeiRz), the equation (4.62) can be written as,
P1 (y(1)) = 11 (4.64)

67

Applying the Laplace transform with respect to one time independent variable yields

and zero initial conditions yields,

PMUYU) = U1 (4.65)

where Y“) and U0) are the Laplace transform respect tl of the vectors y( 1) and u

and the dynamic matrix

2
251 + 251 +12 —2s, — 4

 

 

 

 

P = 2
W) —231 — 4 51 + 231 + 4
Solving for Y“) yields,
Y0) = G1,(1)U1 (4.66)
where
v- 2 .1
G _ 2511+6szll+203f+16s1 +32 2s?+6s¥+203¥+1661+32
1’0) 2314-4 2314-231 +12
4 3 2 4 3 2
L 231+631+2031+1631+32 231 +681+2081+1631+32 J

 

 

is the linearized matrix transfer function and is equal to the inverse of the linearized

dynamic matrix P1

(1)'
The term G2 (U2), is found by equating the elements multiplied by the coefficient
a2 in (4.61),

2 0 -- 2 —2 . 12 4
lo llwl—2 214441-.. .41,» (4.6-5
Using the functional (463), equation (4.67) can be written in the form,
P1 (Y(2)) = 0 (4.68)

Taking the laplace transform of (4.68) for zero initial conditions yields,

PIMYQ) = 0 (4.69)

68

 

For nonsingular P1.(2)’ Y(Q) = 0. This TGSUIt, Y(2)

nonlinear system is an odd function. This result shows the second Volterra operator,

2 O, is expected because the

G2(U2) = Y = o (4.70)

The term G3(U3), is found by equating the elements multiplied by the coefficient

a3 in (4.61),
20.. 2—2, 12-4 80 3
[0 1])f3) + {—2 2 Jy(3)+ [_4 4 :lyE3) + l: O 0 ] Dbl—O (4.71)
Using the functional (4.63) and substituting (4.50) into (4.71) yields,

8 0 -
p1 (y(3)) = - l 0 0 ] y(3) (4'72)

Applying the Laplace transform to (4.72) with respect to three time independent

variable yields,

P Y

8 0 ~
1,(3) (3) = - l ] Y (4-73)

0 O (3)
Substituting the multiplicative operator (4.54) into (4.73),

Pl’(3)Y(3) = —P3¢ (C(1)(U1)a G1(U1)2G1(U1)) (474)

8 0

h P:
were3 [00

J and

2
P = [ 23(3) + 25(3) +12 —23(3) — 4
15(3)

2
—-2s(3) —4 5(3) +2s(3)

and substituting (4.66) in (4.74) yields,

+4

 

Solving for Y by inverting P1

(3) (3)

Y = —G1,(3)P3¢ (G1,(1)U,,G

(3) U1, G1,1U1) (4.75)

1,(1)

where G143) = [P1,(3)]—1

The Volterra model truncated at the third expansion includes equations (4.66) and

. Equation (4.75) is the third term in the Volterra model.

(4.75) (Figure 4.6),
Y1—3 = G1(U1) + G3(U3) (4-75)

69

 

 

 

 

 

 

 

 

 

 

 

u y
‘ Linearized -€\ '
Model , _ ‘
u
2 GI (O) yz
—EE-’ P3 . Linearized
G3 (°)

 

 

 

Figure 4.6. Simulink Diagram for the Volterra Model Including the First and the Third

Expansions

The term G 4(U 4), is found by equating the elements multiplied by a4 in (4.61),

20..+2—2,+
013(4) 423(4)

[12 —4]y +[8 0“! 3(y(1).1):y(2)11 i=0
_4 4 (4) 0 0 13(y(1).2) y(m J

Using the operator (4.64), equation (4.77) can be written in the form,

.4458 mm)

0 0 (3(3’012) y(2).2 1

Using (4.70)
11
(2)71 _
Y(2) [ y l — 0
(2).?

P1 (Y(4)) = 0 => G4(U4) = 0

into (4.78),

Result (4.80) is expected because the nonlinearity is an odd function.

(4.77)

(4.78)

(4.79)

(4.80)

The term G5 (U5), is found by equating the elements multiplied by the coefficient

70

a5 in (4.61),

2 0 2 2 12 4 '(3 (y ):y i
[0 1]y(5)+[_2 2ly(5)+l— 4 4l(5) “Lil 0:“ (1)1 (mlzo (4'81)
1 3,2(y(1)2) y(3).21
Defining the vector,
2
r l 1
~ = y(l),ly(1),1y(3).1 1 = (you) 9(3),1 482
y(21011) ly y y i l 2 l (’ )
(1)12 (1)12 (3)2 ( (yaw) 47(3),, ,

Substituting (4.82) into (4.81) and using the operator in (4.63) yields,

24 0 ..
p1 (315)) = — l 0 0 ] 3(2,0,1) (4'83)

Applying the Laplace Transform to (4.83) respect five independent time variables

yields,
P11(5)Y(5) = _3P3Y(21011) (4'84)
where 2
P _ 28(5) + 23(5) +12 —2s(5) — 4 (4 85)
11(5)_ —23 —4 52 +23 +4 '
(5) (5) (5)

Substituting the vector

 

 

.. Y Y Y
1 z (111) (11) (111) -111 ,. ,1 . 11111
(210.1) (y ) (y ) (Y ) (1) (1) (3 )
(1)12 (1)12 (3)12
and inverting the matrix P1 (5) into (4.84) yields,
Y(s) = )P(15¢(Y )’(Y(1)’Y 3)) (4.87)
—1
where G1,“) = [P1,(5)l and P5 = 3P3 and previous results (4.66), and (4.75) yield,
Yo): G111(1)U
Y(3) = ‘G11(3)P3¢ (G1 (1)U1 G1 (11)U1 G1 1(11)U )

71

all terms of the fifth order expansion can now be computed from the linearized transfer

function G110!)

The Volterra model truncated at the fifth expansion, includes equations (4.67),

(4.77) and (4.87) (Figure 4.7)

 

 

 

 

 

 

Y1-5 = G1(U1) + G3(U3) + G15(U5) (4.88)
u Volterra Model ’ y
1 ; Including 1't and 3"1| > . 1 >
5 Expansions ‘%
u2 (Figure 4.6)

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

(0'50)

 

 

 

Figure 4.7. Simulink Diagram for the Volterra Model Including the First, the Third and
the Fifth Expansions

The time responses of the truncated Volterra models in (4.66), (4.75) and (4.87)
are compared to the response of the nonlinear ODE model (4.60). The two inputs
used to excite the system are sine waves with amplitudes -0.3 and -0.15 at frequencies
of 2.27 rad/sec and 1.5 rad/sec respectively. These frequencies are the two natural
frequencies of the linearized system and are determined by calculating the roots of the

characteristic equation of G in (4.66). The linearized model exhibits the greatest

11(1)
error between all the three responses. The error shown by the truncated Volterra

72

model that includes the linearized and the third expansion terms is substantially
smaller than the error exhibit by the linearized system (Figures 4.8 and 4.9). The
error exhibit by the truncated Volterra model that includes the linearized, the third
and the fifth order expansion term is the smallest of all the three responses. For
a specific bound at the input, the more expansion the Volterra model includes, the
smaller is its error. Theoretically, if the Volterra series is convergent, the error between
the Volterra model and the nonlinear ODE will tend to zero when an infinite number

of expansions is used .

4.7 Summary

In this chapter a method to generate frequency domain Volterra system models from
MIMO external port-based ODEs is presented. Even though a similar procedure
for the SISO case is available in the literature [10], a method for MIMO port-based
external, nonlinear systems has not been published. The MIMO procedure requires
a nonlinear vector operator to obtain the frequency domain kernels. This nonlinear
operator is not required in the SISO case. An advantage of using this MIMO method
is that the two standard MMM nonlinear formats can be easily obtained. The first
format is the Volterra transfer function representation. This format is traditionally
used in model analysis and simulation. The second format is the Volterra dynamic
representation. This format is used in the assembly of nonlinear system models. The
Volterra dynamic representation was recognized and defined for first time in this work.
The Volterra dynamic representation is important in the MMM because this external
port-based model format protects internal design details. This procedure generates

the Volterra dynamic port—based system models in an explicit form.

73

 

Y(llflfil’dfi). ODE

WHO)

 

Figure 4.8. Nonlinear ODE response y, (t) compared with different truncated Volterra
transfer function responses. y(l) is the linearized response, y(3) is the third expansion and
y(5) is the fifth expansion.

74

Q

 

Figure 4.9. Nonlinear ODE response y2 (t) compared with different truncated Volterra
transfer fimction responses. y(l) is the linearized response, y(3) is the third expansion and
y(5) is the fifth expansion.

75

CHAPTER 5

Assembly of Nonlinear Physical Systems

5.1 The MMM Algorithm to Assemble Nonlinear Physical
Models

The algorithm for assembling nonlinear dynamic models through a networked envi-
ronment (Figure 2.1), is explained in this section. Two types of information and a
procedure are required to assemble nonlinear physical models. The first type of infor-
mation is the component connectivity structure of the assembly. The second type of
information is the Volterra dynamic model for each of the components that constitute
the assembly. Once these two sets of information are available at the assembly agent,
the standard procedure to assemble the component models is initiated.

The component connectivity structure of an assembly is enclosed in its constraint
matrix S. The matrix S is dependent on the order of connection between components.
This matrix equates each component output in the unconstrained output vector Yc
to an assembly output in the assembly vector Ya based on the physical constrained
associated with the assembly. Matrix S should be available in the assembly agent
before the assemble procedure is initiated .

If the 2th component model includes q expansion terms, the component Volterra

dynamic model has the form,

{P2,1,(v)1Pi,2,(v)1' " 1Pi,q,(v)} (5-1)

76

For any dynamic matrix Pi, j1(v)’ the index 2', j and 1) respectively indicates the compo—
nent number, the expansion term and the number of independent complex variables

The standard procedure to assembly Volterra dynamic models include 2 steps. In
the first step, the Assembly Agent uses all the Volterra dynamic component dynamic

models to build the unconstrained Volterra dynamic model. This model has the form,

{Pc,1,(v)1 Pc,2,('v)1 ° ' ' 1 Pc,q,(v)} (52)

where (v) is the number of independent complex variables. In an assembly constituted

by k components with a total of r ports, \7’(1 g i S q), each (1‘ x r) unconstrained

component matrix in (5.3) has the form,

F P111102) 0 O 0 q
0 P2,,-,(,,) 0 0

0 0 '°. 0
1. O 0 0 Pk,’i,(’U) .(

Pc,z',(v) = (.33)

 

 

In the second step, the assembly agent generates the Assembly Volterra Dynamic

model. This model is,

{Pa,1,(v)1Pa,2,(v)1' ' ° 1 Pa,q,(v)} (5'4)

Each assembly Volterra Dynamic Matrix in (5.4) is determined using the transforma—
tion,

Pa,i,(v) = STPc,i,(v)S (5-5)
The reason of using (5.5) to determine the Assembly model is explained in detail at
the Appendix B.

The assembly model (5.4) is in the same format than the component model (5.1)
and can be recursively used as a component model to built higher order physical
assembly models. The Volterra dynamic model representation is particularly effective
as the standard assembly nonlinear MMM format because constraints (2.19a) and

(2.19b) can be easily applied [38]. In contrast, the Volterra transfer function model

77

(4.15) cannot be used as the assembly nonlinear MMM format because (2.19a) and
(2.19b) cannot be substituted directly [38]. Even though the Volterra dynamic model
format is effective for the nonlinear model assembly process, it is not appropriate for
simulations because is not possible to use it to solve directly for the outputs given

the inputs. A different model format is required for simulation.

5.2 Simulation of the Assembled Nonlinear Physical Model

The Volterra transfer function model representation (4.15) is chosen as the nonlinear
MMM simulation format because it can be directly used to solve for the outputs given
the inputs. The simulation of the assembled nonlinear physical model is done in two
steps:

In the first step, the Dynamic Assembly Matrix P011101) is inverted to generate
the Volterra Transfer Function Matrix of the assembly GOA“).

In the second step, the assembly Volterra Transfer Function Matrix Ga,1,(v) and
each of the Volterra Dynamic Assembly Matrices P03“), - - - 1Pa,q1(v) are substituted
with the appropriate index (2)) into the standard Volterra Transfer Function formats

(4.17a), (4.17b), (4.17c),- - - ). For the first three standard formats,

Ya,(1) = Ga,1,(1)U1 (5-63)
Ya,(2) = ‘Ga,1,(2)Pa,2,(2)¢ (Ya,(1)1Ya,(1)) (5-6b)

Y411(3) = —Ga111(3) l Pa131(3)¢ (Ya1(1)1Ya1(1)1Ya1(1)) + Pa,21(3)¢(Ya,(1)1Ya1(2)l l
(5.6c)

The Volterra transfer function of the assembly is,
Ya = Ya,(1) + Ya,(2) + Ya,(3) + ' ' ' (5.7)

A simulation of a Volterra model truncated up to the second expansion [10] uses

recursively equations (5.6a ) and (5.6b). Initially the matrix Ga,1,(1) is multiplied by

78

the input vector U(sl). The resultant output is the vector Ya,(1) in function of the
complex variable 31. The inverse Laplace transform respect 31 is applied to Ya,(1)
what yields the the degree-1 time homogeneous output ya,1(t1). The time variable
t1 is redefined as t = t1 what finally yields the degree-1 time homogeneous output
ya,1(tl-

Equation (5.6b) operates the previously calculated vector Ya“) into
¢(Ya,(1)1 Ya,(1)). The resultant vector output of this operation is multiplied by the
matrix P0342), and then by the matrix —Ga,1,(2) what yields the vector Ya,(2) which
is a function of the complex variables 31 and 32. The inverse Laplace transform is
applied sequentially to Ya,(2) two times, first respect 31 and then respect to 32 what
yields the the degree-2 time homogeneous output ya,1(t1,t2). The time variable t1
and t2 are redefined as t = t1 = t2 what finally yields the degree—2 time homogeneous
output ya,2(t).

The time response of the Volterra model truncated up to the second expansion
ya,1-2(t) is obtained by adding the degree-1 time homogeneous output to the the

degree-2 time homogeneous output what yields,

Ya1_2 (t) : Ya,1(t) + Ya,2 (t) (58)

Simulation of higher order terms are obtain following a similar procedure.

5.3 Example

The assembly of two component models using the nonlinear MMM algorithm is de-
scribed in this section. The first component is a linear model with one external port.
This port is the input-output pair (143,3;3). The second component is a nonlinear
model with two external ports. These ports are the input-output pairs (ul, yl) and
(142,312). Two joins are used in the assembly. The first joint connectes the port (213,373)

of the linear component and the port (212,312) of the nonlinear component to an as-

79

sembly port ((12, 372). The second join connects the component port (111,311) to the
assembly port (211,311). The assembly model has two ports, these are the input-output
Pairs (1111171) and (1721172) (Figure 5-1)-

J71 I a) 5" 2 luz

‘ “1 Nonlinear “2 ‘ “3 Linear
y Component y2 y3 Component

1

 

 

 

 

 

 

 

 

 

 

 

 

Figure 5.1. Port Connection Diagram of the Assembly

The linear model describes a one-port linear mass—spring system (Figure 5.2), with
mass m3 and spring constant 163. The port uses the input—output pair (113,313) where
(43 is the input force applied to the mass m3 and 313 is its output displacement. The

ordinary differential equation that represents this linear mass-spring—system is,

m3g'j3 + k3y3 = 213 (5.9)

 

 

Y3“)

Figure 5.2. Linear Mass-Spring System

The nonlinear model describes a two-port mass—spring-damper system (Figure
5.3) composed of two masses, one nonlinear spring, one linear spring and one linear
damper. The mass m1 is connected to a fixed point through the nonlinear spring with
force fN L = k1(y1 + yi). The mass m2 is connected to the mass ml with a linear
damper having damping coefficient ()1 and a linear spring with constant k2. The first

port uses the input-output pair (111,371) where u1 is the input force applied to the

80

mass m1 and y1 is its output displacement. The second port uses the input-output
pair (ug, y2) where 142 is the input force applied to the mass m2 and y2 is its output
displacement.
1110) “20)
-—'> k —>
2

Ngglrlfilegar _ ‘flllllllllllllllll— m2
b1 L.>

y1(t) 142(1)

 

 

 

 

 

 

 

 

 

 

 

 

Figure 5.3. Double-Mass Spring Damper—Nonlinear System

The ODE that represents this nonlinear mass-spring-damper system is,

mlyl + bl(yl ‘ 392) + klyl +1913/13 + k2(yl _ 112) = "1
T”2132 + b1(y2 " 91) + k2(y2 ‘ 311) = “2

(5.10)
The assembly is performed by attaching the mass m3 of the linear system to the
mass m2 of the nonlinear system (Figure 5.4). The resultant assembled system has two
ports. The first port uses the input—output pair ((11, m) = (ul, yl). The second port
uses the pair (112, 172) where 112 is the input force applied to the mass m0 = mg + m3
and 372 is its output displacement. Traditionally and ad-hoc procedure is used to
obtained the ordinary differential equation of the assembly. This procedure is difficult
to automatize because it usually requires a new reformulation using physical laws or
energy methods. The new ODE for the assembled system is,
miljl + b1(?21—§2)+ 1311.714" 1911713 + k2(371“ 92) = {‘1 (511)
(m2 + m3)g2 + (310.72 _ .01) + k2<172 _ 171) + [93.732 = {‘2
The nonlinear MMM algorithm requires the Volterra dynamic model of each as—
sembly component. A standard procedure is used to derive the corresponding Volterra

dynamic models of the linear and nonlinear ODEs in (5.9) and (5.10). Each Volterra

dynamic model provided has a finite number of expansion terms.

81

u] -4. 1102 = 112 + U3

     

Figure 5.4. Assembly of a Linear and a Nonlinear Mass Spring System

The third order truncated Volterra dynamic model for the nonlinear ODE in (5.10)

has the form,

{P1,1,(U)’P1,2,(‘U)’P1,3,(‘U)} (5-12)
where
2
1 1( ) = [ "mm + blsm + k1 + k2 21,18”) _ k2 ] (513a)
, , ’U
_bls(v) — k2 mzsw) + (218(1)) + kg

0 0
P112,(v) = [ 0 0 ] (5.13b)

161 0
P134”) = [ 0 0 :l (5.13C)

The third order truncated Volterra dynamic model for the linear ODE in (5.10)

has the form
2

P211102) = m3s(v) + k3 (5.14a)
P 2,2.(11) = P2,3,(v) = O (5.14b)

The matrix S is generated by relating the component and assembly variables. The
component variables are y1(t), y2(t), y3(t), u1(t), u2(t) and 113(t). The assembly vari-
ables are 391 (t), 372(1)), {q(t) and 112(t). The two constraint equations for this example

are:

3,; (t) 1 0 -
1120) = 0 1 [151(1)] (5.15a)
.7130) 0 1 11205)

1 O {u(t) U10)
0 1 [ ]: U2“) (5.15b)
0 1 “2m 1130*)
where,
1 0 13(5) U1(s) ~ “
S: 0 1 ,Yc(s)= 14(3) 11%“): U2“) ’Y“(S)=lf1l:fl’U“(S)=lglf:il
0 1 Y3(8) U3(3) 2 2

In the first step the Assembly agent generate the unconstrained Volterra dynamic

model. This model have five matrices which are,

I l
i , [ P1,2,(v) 0 ] , [ P1,3,(v) 0 J} (5.16)

0 0 0 0
In the second step, the Assembly Agent determines the five dynamic assembly

P111101) 0
0 P211411)

 

 

matrices. These are:

 

 

 

P O
_ T 1111(1))
Pa’1’(v) _S 0 132.1 (v) S
2 ’ (5.17a)
P _ m18(v) + b18(v) + [61+ 162 —b18(v) — k2
a111(v) —b18(v) — k2 (m2 + m3)s%v) + (918(1)) + k2 + [133
P 0 0 0
Pa,2,(v) = ST |: 131(1)) 0 ] = [ 0 O :l (517b)
P O k 0
T 1,3, _ 1
Pa,3a(U) = S l: 0 (U) 0 J S _ [ O 0 :l (517C)
Finally the Volterra dynamic model is,
{Pa,1,(v)1Pa,2,(v)1Pa,3,(v)} (5-18)

The simulation is done in two parts

1. The dynamic matrix P011100 is inverted to generate the assembly transfer function

 

 

matrix,
(7712+7n3)s%v)+b1S(,,)+k2+k3 b1 5(1’)+k2
G _ [P ]—1 _ as4+bs3+032+ds+e 334+bs3+csz+ds+e
“114”) _ “114”) _ b k m 52 +b s H: +k
18(1))+ 2 104 1(6) 1 2
asJ+bs3+csg+ds+e 334+bs3+csz+ds+e

 

 

(5.19)

83

 

where

a = m1(m2 + 7713)
b = b1(m1 + mg + m3)
C = (k2 + k3)m1 + (k1 + k2)m2 + (k1 + k2)m3
d = b1(k1 + (£3)
6 = [£11162 + (131,63 + k2k3

2. The assembly transfer function matrix (30,130,) and the four dynamic assembly ma-
trices Pug”), - -- ,Pa,5,(v) are substituted into the standard Volterra transfer func-

tion formats. The first five Volterra formats for this system are,

Y(l) = (311,1,(1)U1 (5.20a)

0 0
Y(2) = l O 0 J (5.20b)

0 0
Y(4) = l 0 0 ] (5.20d)
Y(5) = 3Ga111(5)Pa13(5 5) w(Y(1)1Y(1)1Y(3)) (5.208)

The truncated Volterra model up to the fifth expansion is,

Y1—5 = Y(1) + Y(3) + Y(5) (5.21)

One way to simulate the Volterra transfer function model is using a Simulink dia-

gram. The Simulink diagram for the truncated model in (5.21) is shown in Figure 5.5.

The time responses of the linear, third and fifth truncated Volterra models are
compared to the nonlinear ODE of the assembly in (5.11) using the system parameters,
m1 = 2, m2 = 1, m3 = 1, bl = 2,131: 8, k2 = 4 and k3 =1. The two inputs used
to excite the system are sine waves with amplitudes -0.3 and -0.15 at frequencies
of 2.36 rad/ sec and 1.31 rad/ sec respectively. These frequencies are the two natural

frequencies of the linearized system of the assembly and are determined by calculating

84

 

 

 

 

 

 

 

 

 

 

 

 

Ga,l,(3)

 

 

 

 

 

 

 

 

 

 

 

 

Pa,3.(3)

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

. G410)

 

 

 

 

 

 

 

 

 

 

___1 -

Ga,l,(3)

 

 

 

 

 

 

 

 

 

V

 

 

 

 

 

 

 

 

 

 

 

 

 

 

IE]

 

 

 

Input
Input

:3“

Figure 5.5. Simulink Program Used to Simulate the Volterra Model

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

the roots of the characteristic equation of G110) in (5.19). The linearized model
exhibits the greatest error between all the three responses. The error shown by the
truncated Volterra model that includes the linearized and the third expansion terms is
substantially smaller than the error exhibit by the linearized system (Figures 5.6 and
5.7). The error exhibit by the truncated Volterra model that includes the linearized,
the third and the fifth order expansion term is the smallest of all the three responses.
For a specific bound at the input, the more expansion the Volterra model includes,
the smaller is its error. Theoretically the error between the Volterra model and the

nonlinear ODE will tend to zero when an infinite number of expansions is used.

86

y(1)+y@)
\

 

Figure 5.6. Nonlinear ODE assembled response y, (t) compared with different tnmcated
Volterra transfer function responses. y(l) is the linearized response, y(3) is the third ex-
pansion and y(5) is the fifth expansion.

87

 

 

Figure 5.7. Nonlinear ODE assembled response y,(t) compared with different truncated
Volterra transfer flmction responses. y(l) is the linearized response, y(3) is the third ex-
pansion and y(5) is the fifth expansion.

88

5.4 Summary

In this section an original method to assemble physical MIMO physical external port-
based Volterra system models is presented for first time. Volterra representations are
appropriate to the MMM approach to distributed model assembly because they pro-
tect internal design details. This method extends the MMM approach to nonlinear
port-based models through a method that is computationally efficient because it is
recursive and does not require matrix inverses. These properties make this Volterra-
based assembly procedure ideal for distribution and assembly of nonlinear physical
engineering system models. This method does not present a new modeling methodol—
ogy but rather a new strategy to distribute and assemble existent nonlinear models.

The main purpose of this nonlinear model assembly methodology is to enhance
Global Engineering Design by providing a procedure that satisfies the four require-
ments of a Global design environment. These four requirements are: standardize
model components, single query exchange of model information, external input-output
model formats and a recursive assembly.

An example of an assembly of one, nonlinear 2-port mass-spring—damper system
and one, linear 1-port mass-spring system was presented. The resultant Volterra
model of the assembly includes the first five expansions. Different truncated MIMO
Volterra models for the assembly where also simulated using Simulink, and their
responses were compared with the nonlinear ODE of the assembled system. The
response results show increasing simulation accuracy as the number of expansions
included in the Volterra model is increased. Experience with these simulations models
shows convergence a requires bound on input amplitudes. Convergence radii are
available for SISO Volterra models [13]. Future work will be needed to find these

limits for the MIMO models derived from the MMM assembly methodologr.

89

CHAPTER 6

Conclusions

The nonlinear Modular Model distribution and assembly methodology presented in
this work, provides six (6) significant, original contributions to mechanical engineer-
ing. These contributions are:
1. An extension of the Modular Modeling Methodology (MMM) to nonlinear systems.
2. The recognition that two model formats are necessary for physical system model
assembly and simulation.
3. A method to assemble physical port-based affine ODEs around an equilibrium
operating point.
4. A method to obtain subsystems operating points from the system assembly oper-
ating point outputs.
5. A method to obtain frequency domain Volterra system models from MIMO port-
based ODEs.
6. A method to assemble physical MIMO port-based Volterra system models. Each
of these contributions is individually important to the development of a Global Engi-
neering Design (GED) strategy.

This methodology has two (2) limitations. These limitations are:
1. The nonlinear subsystems must be modeled using port-based Volterra representa-
tions.

2. The assembled system model is valid if its subsystem models are valid. These

90

limitation still allows the model assembly of a high variety of nonlinear port-based
physical system models.

The nonlinear model distribution and assembly methodology developed in this
work enhances Global Engineering Design by providing a procedure that satisfies the
four requirements of a Global design environment: a standard model format, single
query exchange of model information, external input-output models and a recursive

assembly method.

91

 

APPENDIX A

Singular Internal Stiffness Matrix

A singular internal dynamic matrix Piz-(s) occurs when the assembled system includes
internal disconnected subsystems or non-observable modes. When Pii(s) is singular,
a modification of the approach taken in (2.44) it is required. The objective is to
remove the internal singular degrees of freedom associated with the disconnected
subsystems. The internal singular degrees of freedom are removed using a null-space
transformation. Once these degrees of freedom are removed, the modified internal
dynamic matrix Pii(s) is non-singular and the condensation continues.

Initially the rank p orthonormal null space T(s) (v x p) matrix of Pz-i(s) is deter-
mined. The columns of this matrix are the p eigenvectors of Piz-(s) whose eigenvalues

are zero. The null space T(s) defines the (v x v) transformation matrix

e(s) = [ (I, he] (A1)

constructed with the (2) X v) null space T(s) and an arbitrary set of independent
vectors. The independent set vectors constructed here are an ((2 — p) x (v — p)
identity matrix filled below with a (p x p) zero matrix.
The above null space transformation is used to define a new set of internal output
variables
Y1(8) = Q(8)Y(S) (A?)

where the new internal output variables Y(s) = [ 23(9) Y0(s) ]T are the (v— p) non-

92

Singular output 21(5) combined with the p singular outputs Y0(s). Work must be
conserved under this coordinate transformation. Defining the transformed generalized

input U(s) and applying a work analysis parallel to (2.16) above,

Y.(s)TU1-(s) = Y(s)TU(s) (11.3)
Substituting (A.2) into (A3) yields
Q(s)TU1-(s) = 11(3) (4.4)

Applying this transformation symmetrically to (2-5-11) yields

P6603) Pei(8)Q(S) ] —%:%* =[ Ue(5) ] (A.5)

T8 '8 T8 "3 8 T8 '3
Q ()P..() Q ()P..( )Q() 170(1) Q ()U.()

Transformation Q(s) contains the null space transformation with Pii(s)T(s) so that

(A.3) becomes

366(3) 1:362 (5) 0 Y6“) [36(3)
Pze(8) P 24(5) 0 Yi(5) = U4(8) (A-G)
0 0 Y0(8) U0(8)

None of the system equations involve Y0(s) and this vector can be removed from
the system of equations. The last p rows and columns of the transformed system are
now removed and the analysis continued with LIZ-(s) = 0 applied on the transformed

matrices to find

Pe(s)Ye(s) = Ue(s) (14.7)

where

93

APPENDIX B

The Assembly Operator

The reason of using (5.5) to obtain all the assembly Volterra dynamic matrices is
explained here. Since the matrix Ga1l1(v) is the linearized transfer function of the

assembly, its inverse, Pa,“

.0), is the linearized dynamic matrix of the assembly. It

was shown previously that Pa,1,(,,) can be obtained using,
T
Pa1l,(v) = S Pc,1,(v)S (B.1)

The dynamic matrix Pa,1,(v) satisfies (4.20). Applying the operator Dy,y,y1'" (e) to
both sides of (B.1) yields,

Pa121(v) 2 BMW (P1111410) = Dyssa'" (STPc111(v)S) (B2)
The matrix of constants S can be taken out of the operator to yield,
1102),) = STDyW,.. (Pc,1,(,)) s (3.3)
Using (4.20) in the left side of (B3) yields
Pa,2,(v) = STPc,2,(v)S (13-4)

Applying again the operator D .. (o) to both sides of (8.4) yields

yak?“
__ _ T
Pa,3,(v) — Dy,y,y1'" (Pa,2,(v)) — Dy,y,y1"° (S Pc,2,(v)S) (B5)

The matrix of constants S, can be again taken out of the derivative operator to yield,
_ T
Pump) _ 3 pm,” (1362),») 3 (B6)

94

Using again (4.20) in the left side of (A6) yields
T
Pa,3,(v) = S P431003 (13-7)

Applying sequentially a similar procedure for any i = 3, 4, 5, - - - yields the assembly

operator (5.5) ,

95

BIBLIOGRAPHY

[1] Elmqvist, H., 1999, ”Modelica: A language for Physical System Modeling, Visu-
alization and Interaction”, CACSD’99 IEEE Symposium on Computer Control
System Design, Hawaii, August 22-27.

[2] Reichenbach, D., Radcliffe, C., and Sticklen, J., 2004, ”Modular Distributed
Models of Structural Dynamics”, Automated Modeling 2004 ASME IMECE,
Anaheim, CA., November 13-19.

[3] Tian, G., Yin, G., and Taylor, D., 2002, ”Internet-Based Manufacturing: A re-
view and a new Infrastructure for distributed Intelligent Manufacturing”, Journal
of Intelligent Manufacturing, Vol. 13, pp. 323-333.

[4] Cu, B., and Asada, H., 2004, ”Co-Simulation of Algebraic Coupled Dynamic Sys-
tems Without Disclosure of Proprietary Subsystem Models”, Journal of Dynamic
Systems Measurement and Control, Vol. 126, pp.1-13.

[5] Zienkiewicz, O., 1989, The Finite Element Method, McGraw-Hill, London.

[6] Mattsson, E., 1993, ”Index Reduction in Differential Algebraic Equations Using
Dummy Derivatives”, SIAM Journal on Scientific Computing, Vol. 14(3), pp.
677—692.

[7] Karnopp, Dean C., Margolis, D., and Rosenberg, R., 2000, System Dynamics:
Modeling and Simulation of Mechatronic Systems, John Wiley Sons, Inc., New
York.

[8] Byam, B., 1999, Modular Modeling of Engineering Systems Using Fixed I-O
Structure , PhD. Dissertation, Michigan State University, East Lansing, MI.

[9] Byam, B., and Radcliffe, C., 2000, ”Direct-Insertion Realization of Linear Modu-
lar Models of Engineering Systems Using Fixed Input-Output Structure”, ASME,
Proceedings of DET2000: 26th Design Automation Conference, Baltimore, MD,
September 10—13.

[10] Rugh, W., 1981, Nonlinear System Theory, The Johns Hopkins University Press.

96

[11] Brilliant M, March 1958, ”Theory of the Analysis of Nonlinear Systems”, Tech-
nical report 345, Massachussetts Institute of Technology, Research Laboratory
of Electronics.

[12] Sandberg W. I, June 1983, ”Series Expansions for Nonlinear Systems” Circuits,
Systems, and Signal Processing, vol. 2, no. 1, pp. 77-87.

[13] Boyd, S. and Chua L, November 1985, ”Fading Memory and The problem of
Approximating Nonlinear Operators with Volterra Series”, IEEE Transactions
on Circutis and Systems, Vol cas-32, No 11.

[14] Wang, T., and Brazil, Thomas, 2000, ”The Estimation of Volterra Transfer
Functions With Applications to RF Power Amplifier Behavior Evaluation for

CDMA Digital Communication, IEEE Microwave Symposium Digest, Vol. 1,
June 2000 11-16, pp 425 - 428.

[15] Draganescu, E., and Ercuta A., 2003, ”Identification of Nonlinearities in Anae-
lastic Polycrystalline Material Using Volterra-Fourier Transform”, Journal of Op-
toelectronics and Advanced Materials, Vol. 5(1), pp. 301-304.“

[16] Li R, and Pileggi T. L., 2003, ”NORM: Compact Model Order Reduction of
Weakly Nonlinear Systems”, DAC 2003, Anaheim, California, June2-6.

[17] Giret, A., and Botti, V., 2004, ”Holons and Agents”, Journal of Intelligent Man-
ufacturing, Vol. 15, pp. 645-659.

[18] Takahashi Y., Rabins, M., and Auslander D., 1970, Control and Dynamic Sys-
tems, Addison-Wesley Publishing Company.

[19] Morari, M., and Zafiriou, E., 1989, Robust Process Control, Prentice Hall, En-
glewood Cliffs, New Jersey

[20] Genta, G., 1999, Vibration of Structures and Machines, Springer-Verlag, New
York.

[21] Skogestad, S., and Postlethwaite, I., 2000, Multivariable Feedback Control, John
Wiley and Sons, England.

[22] D’Souza, A., 1988, Design of Control Systems, Prentice-Hall, New Jersey.

[23] Buck, R., Willcox, A. 1971, ”Calculus of Several Variables”, Houghton Mifflin
Company, Boston

[24] Favier, G., Kibangou, A., and Y., Khouaja, 2004, ”Nonlinear System Modelling
by Means of Volterra Models. Aproaches for Parametric Complexity Reduction”,
Simposium Techniques Avancees et Strategies Innovantes en Modelisation et
Commandes Robustes des Processus Industriels, Martigues, September 21-22.

97

 

[25] Bikerlund, Y., and Hanssen, A., 2003, ”On the Estimation of Nonlinear Volterra
Models in Offshore Engineering”, International Journal of Offshore and Polar
Engineering, Vol. 13(1).

[26] Boaghe, O., and Billings S., 2003, ”Sub-harmonic Oscillation Modeling and
MISO Volterra Series”, IEEE Transactions on Circuits and Systems, Vol. 50(7),
pp 877-884.

[27] Storrs, Samuel M., 1999, ”Noise Chacterization of devices for Optical Comput-
ing”, Ph.D. thesis, Texast Techn. University.

[28] Schetzen, M., 1980, The Volterra and Wiener Theories of Nonlinear Systems,
John Wiley

[29] Yangwang F ., and Jiao, L, 2001, ”MIMO Volterra Filter Equalization Using
Pth-Order Inverse Approach” Proc. of the IEEE International Conference on
Acoustics, Speech and Signal Processing, Vol. 1, pp 177-180.

[30] Dobrowiecki, T., P., and Schoukens, J ., 2004, ”Linear Approximation of Weakly
Nonlinear MIMO Systems”, Proc. IMTC 2004, Como, Italy, May 18-20.

[31] Yangwang F., and Jiao, L, Ruixuan, W., Pan. J. 2000, ”Identification of MIMO
Nonlinear Systems Based on Volterra Kernel Matrices” Proceedings. of the 3rd
World Congress on Intelligent Control and Automation, June 28-July 2, Hefei,

PR. China.
[32] Greenwood, D, 2003, Advanced Dynamics, Cambridge University Press.

[33] Nam,C., 2006, Aeroelastic Analysis Using Matlab, at
http://ctaspoly.asu.edu/chnam/ASE Book/

[34] Motato, E., Radcliffe, C. and Reichenbach, D., 2006, ”A Procedure for Obtaining
Multi Input Multi Output Volterra Models From Port-Based Nonlinear Differen-

tial Equations”. Submitted to the 2007 American Control Conference, New York,
NY.

[35] Chua, L., and Lin, Pen-Min , 1975, Computer-Aided Analysis of Electronic Cir-
cuits, Prentice-Hal, Inc., Engewood Cliffs, New Jersey

[36] Iserman, R., 1992, Adaptive Control Systems, Prentice Hall.

[37] Kerr, Brad., 2000, ”Redesigning work Processes and Computing Environments”,
Automotive Engineering International, Vol. 108(7), pp 147-149.

98

 

[38] Motato, E., Radcliffe, C. and Reichenbach, D., 2006, ”Networked Assembly of
Linear Physical Models”. Submitted to the ASME Journal of Dynamic Systems,
Measurement , and Control.

[39] Frank P, and McFEE, June 1963, ”Determining Input-Output Relationship of
Nonlinear Systems by Inversions” IEEE Transactions on Circuit Theory, pp 169—
180.

[40] Christensen G, December 1969, ”On the convergence of Volterra Series” IEEE
transactions on Automatic Control, 13, 736-737.

[41] Czarniak A and Kudrewicz, August 1984, ”The Convergence of Nonlinear Series
for Nonlinear Networks” IEEE transactions on Circuit and Systems, 31, 751-752.

[42] Chatterjee A., and Vyas N., February 2000, ”Convergence analysis of Volterra
Series Response of Nonlinear Systems Subjected to Harmonic Exitation”, Journal
of Sound and Vibrations. 236(2), pp 339—358

[43] Tomilson G. R. and Manson G., February 1996, ”A simple Criterion for Establish-
ing an Upper Limit to the Harmonic Excitation Level of the Duffing Oscillator
Using the Volterra Series”, Journal of Sound and Vibrations. 190(2), pp 751-762

99

 

   

1lull)((1)