Equation-Free System-Level Dynamic
Modeling and Analysis in Energy Processing
Abstract—The paper describes a method for system -level
modeling and analysis in energy processing that does not rely
on equations to capture global dynamics. Instead, it acts directly
on component-level (micro) models, and thus avoids additional
modeling assumptions that are often invoked on the system level
to generate an explicit model. The equation-free framework is
capable of directly answering key analytical questions about
system dynamics, such as steady states, stability and bifurcations.
The procedure is illustrated on a switched power electronic dc/dc
converter.
I. INTRODUCTION
Key analytical challenges in energy processing systems
stem
from their intrinsic temporal and spatial multi-scale nature.
Such systems often span more that one physical domain (e.g.,
fuel cells that amalgamate electrical, thermal and chemical
properties). Even a brief review of more than a century of
sustained effort in physics-based modeling of energy processing
components and systems would reveal that much of this
knowledge is, unfortunately, only available in rigid, inflexible
formats such as proprietary and legacy computer code.
The best models available today are often on the component
(microscopic) level, while the tasks facing the analyst are on
a much coarser, system (macroscopic) level. For example, in
power electronic systems switching is a key feature at the
component (e.g., converter) level, but its details are commonly
abstracted at the higher level of electrical, mechanical,
chemical and thermal interactions. The main tool here has
been averaging and its numerous variants. Averaging is based
on assumptions which are reasonable in most cases, but
seldom checked during actual analysis and simulation of large
systems.
For practical utility, we need a framework that
encompasses
both physical (phenomenological) and behavioral (data-driven)
models. There exists a strong preference for physical
models in the energy processing community, justified by the
non linear nature of energy transformations, by the need to
make predictions based on extrapolations and to organize
the prior information, and by a long history of modeling
components for design purposes. Behavioral models are of a
more recent vintage in energy processing, and correspond to
cases when a physical understanding of underlying phenomena
is unavailable, or is too costly to obtain.
A. Algorithmic Models
In motivating the need to broaden the class models used
in energy system studies, we briefly review a classification of
dynamical system models in three groups, as outlined in [1].
The first paradigm, denoted as Newtonian and established in
the seventeenth century, required modeling a system with a
differential equation , and then explicitly solving that equation,
i.e., providing a formula for the solution . While very successful
in examples it tackled (mostly from physics), it also
proved too restrictive, as there is still no general procedure for
finding closed-form solutions when they exist, and there are
elementary functions of practical importance whose integrals
are not elementary.
The second paradigm emerged from the pioneering
work
of Poincare in the nineteenth century. It also begins with
differential equations, but does not require explicit solutions
- instead, it attempts to provide qualitative information about
system behavior directly. This so called qualitative theory of
differential equations brought forward many advances, including
work in power systems (e.g., structure preserving dynamic
models[2]). Foremost on this list are the studies of equilibria,
stability and effects of parametric changes. However, it is
not practical to always require explicit models in terms of
ordinary differential (ODE) and algebraic (DAE) equations
– for example, in cases when models are encapsulated in
proprietary or legacy computer code.
The third paradigm, algorithmic modeling, focuses
on
utilizing models from heterogenous sources, exemplified by
executable computer code, logic automata and tabulated measurement
data. It is prudent to be cautious when fitting models
to data, especially when models are not based on underlying
physics. A strong reminder is provided by the so-called universal
differential equations [3], capable of approaching with
arbitrary precision any continuous behavior with a model that
depends on only five parameters. The sensitivity of a solution
to variation of these five parameters is extreme, of course,
and the model offers no insight into the phenomenon under
study. Note that algorithmic modeling subsumes a large part
of the Poincare program, as it also aims to answer qualitative
questions directly, while combining geometric and algebraic
methods. Models that exist in the form of proprietary or legacy
code naturally belong to this class.
B. Multi-scale Nature of Energy Conversion
One common feature of energy processing systems is their
multi-scale nature. This is manifested in large utility power
systems, where time-scales of interest for control (see Fig. 1,
adapted from [4]) range from intra-cycle (60 Hz) behavior of
power electronic converters, to electromechanical transients
(on the order of seconds) to long-term dynamics involving
power plants (tens of seconds to minutes). Analogous decompositions
are relevant for power electronic converters switching
at tens of kHz or more - from inter-cycle behavior important
for closed-loop control, to studies of behavior of connected
converters that involve phenomena that are several orders
of magnitude slower. Similarly, in industrial electric drives
there is typically a wide sepa ration between dynamics of the
electronic converter and the dynamics of the mechanical and
thermal subsystems.
In conventional, equation-based approaches the transition
from more to less detailed models (model reduction) is
achieved analytically, within the overall modeling task. Such
procedures are based on variants of averaging, and the aim
is to disregard the details (“ripple”) while retaining the parts
most relevant for slower dynamics. Because of these a priori
modeling decisions, the procedure necessarily relies heavily
on assumptions about the regularity of the responses of the
fast modes (e.g., regarding the frequency content of signals
of interest, or balance among phases). These assumptions are
reasonable in a majority of applications, but not necessarily in
all cases, especially when seemingly irregular or anomalous
behaviors lead to cascades of events. Standard modeling
assumptions may fail for several reasons:
• new types of subsystems that violate standard modeling
assumptions (e.g., power electronic converters or distributed
energy sources),
• models of subsystems in nonstandard formats -
measurements-based multidimensional tables or computer
code involving logical variables used for power
system protection - for which it may be hard to apply the
conventional modeling procedures devised for smooth,
differential equation-based dynamics,
• uncertainties in parameter values.
When focusing on standard transient stability models on
the
power system level [4], [5], we observe a number of modeling
idealizations: 1. stator and network transients are neglected,
which is compatible with the presence of only the fundamental
phasors in all components, 2. assumed either balanced operation,
or decoupled positive , negative and zero sequences (for
example, this is only approximately true in electric machines
[6]), 3. decoupled real and reactive loads, 4. simple mechanical
models, including constant power input and rigid body motions
for rotors, 5. simple models for magnetic saturation and for
effects of frequency variations. These assumptions are reasonable
in a conventional framework, as they enable descriptions
of large systems for which explicit ODE or DAE models are
needed. Note, however, that these assumptions may be violated
in many practical cases (e.g., aircraft and ship-based systems
and microgrids).
For these reasons we may be unable to produce an accurate
system-wide (“top-level”) model that would be written in
the standard format of explicit differential-algebraic equations
(DAE) in, say, continuous-time, while at the same time delivering
the needed fidelity. However, we may still be able to
answer some key questions regarding, e.g., stability by acting
directly on a hierarchical multi-scale model. In doing so, we
can use an approximate system-level model as well (when
available), to improve numerics and to speed up the modeling
procedure. On the conceptual level, we are making a transition
from analytical solutions (of explicitly available equations) to
computer–assisted solutions.
II. EQUATION-ASSISTED MODELING
A. Modern Multi-scale Methods
A wider background for our approach is provided by the
heterogeneous multi-scale method [7] and by the “equationfree
modeling of complex systems ” [8], [9], [10]. Note that
traditional system-level (macro-scale) models in power system
have their limitations, most notably accuracy and the empirical
nature of some component models (e.g., loads). At this point
one may be tempted to switch to detailed (micro-scale) models
(EMTP-type in power systems) for system-level studies.
However, this strategy is suboptimal for two reasons: 1) such
models may be unwieldy because of size and complexity (even
with advances in computing technology), and 2) the answers
would likely contain too much data that is of little interest,
complicating the task of extraction of useful information even
further. Thus the basic task of multi-scale modeling is to blend
the two scales in designing computational methods that are
much more efficient than solving the full microscopic model
and at the same time yield system-level information with
the desired accuracy. Our methodology is data-driven in the
sense that the linking between models at different scales is
achieved through the data (in contrast with requiring analytical
solutions).
B. Equation-Free Analysis of Complex Systems
A persistent feature of many complex systems is the
emergence
of system-level coherent behavior from the interactions
of lower-level components (“agents”) among themselves and
with their environment. In standard engineering practice, the
system-level models are often in the form of conservation
laws (energy, momentum, charge or mass) closed through the
constitutive equations of components (e.g., “v-i characteristics”).
In many problems of current modeling practice across
engineering, the models are well developed on the component
level (and often are physics-based), but the closures required
to translate them to accurate system-level descriptions are not
available. This is particularly true in the case of models using
proprietary or legacy computer code. The main idea proposed
in [8], [9], [10] is to consider the component-level simulator
as a computational experiment that one can set up, initialize
and run at will. The hope is that the results so obtained will
allow the analyst to estimate the same information about the
system-level model that would be evaluated from (unavailable)
explicit formulas. This framework extends readily to the case
of hardware-in-the-loop simulations as well.
Short experiments involving the micro model, at states
prescribed by the main, macro level algorithm, give values to
residuals and Jacobians acting on vectors of particular interest;
in a way, traditional numerical algorithms become protocols
for the design of physical and numerical experiments.
In equation-free computation we do not have a subroutine
containing the macro model, but we do have a subroutine
that can evolve the micro model, and we can initialize the
micro model consistently with values of its corresponding
macroscopic variables. When the system-level algorithm asks
for specific Taylor series information of the unavailable macro
model at a particular state point, we then initialize the micro
scale code consistently, perform possibly several short bursts
of computational experimentation with it, and estimate the
numbers of interest. These numbers are then fed back to the
system-level algorithm which can now continue its procedure.
This multi-scale procedure is a in many ways analogous to
matrix-free methods of iterative linear algebra [11]. To perform
tasks such as finding leading eigenvalues or solving linear
equations we do not need the matrix itself. It is enough to
have a sequence of matrix-vector products for a sequence of
known vectors.
We can also use approximate models available at the system
(macro) level to complement computations of the component
(micro) level models. This is the “equation-assisted” framework
proposed for bifurcation analysis of partial differential
equations (PDE) in [12]. The role of the approximate model
is to accelerate the convergence of micro-level calculations
by constructing preconditioners for iterative numerical linear
algebra involved in bifurcation and stability computations.
III. TOOLS FOR ALGORITHMIC MODELING
A. Analysis of Algorithmic Models
Note that the equation-assisted algorithmic modeling
framework
is very flexible. It can be used not only for temporal
evolution of the system-level model, but also for determining
its steady states and other fixed point computations. In particular,
a number of bifurcation studies have been reported in the
equation-free framework [13], [14]. The key for such analyses
is the Recursive ProjectionMethod (RPM) of Shroff and Keller
[15] and its more recent extensions [16]. An extension of
RPM to DAE systems is presented in [17]. In the case of
stability and bifurcation analyses, the key idea is to develop
an efficient single-shooting based technique in which the
explicit calculation of the Jacobian is avoided. The approach
requires only the calculation of the action of the Jacobian
on a low-dimensional subspace corresponding to dominant,
stability-determining eigenvalues and Floquet multipliers . In
geometric terms, matrix-free, timestepper based eigensolvers
are used to estimate the slow eigenvalues and the corresponding
eigenvectors for the timestepper, which should be tangent
to the slow manifold (embodying the missing accurate macro
model). Gaps in this spectrum and the components of the
corresponding eigenvectors are used to locally parametrize the
low-order manifold.
Given the current state of the system x, the time-stepper
computes the future state Φ(x; τ) where τ is the “reporting
time”; note that x(τ) ≡ Φ(x; τ). For later use we introduce
functionΨ(x) ≡ x−Φ(x; τ) whose zeros correspond to steady
states of the original system. In calculating the fixed points of
Ψ(x) one commonly uses iterative algorithms like Newton-
Raphson. The key step in our calculations will be to opt for
matrix-free solvers of linear equations inside such algorithms,
e.g., GMRES, [14]. In doing so, we will avoid the need to
explicitly calculate the Jacobian DΨ ≡ ∂Ψ(x)/∂x. We will be
required to calculate matrix-vector products of this (explicitly
unknown) Jacobian with known vectors v; this quantity can be
estimated using the finite difference approximation DΨ· v ≈
[Ψ(x + εv) −Ψ(x)]/ε for a suitably small ε.
Another key issue is the selection of appropriate systemlevel
states x for which we are trying to extract the algorithhmic
model through bursts of micro-level simulations. We
are helped in this quest by the long history of physics-based
modeling in energy processing systems, which offers valuable
guidance in this respect. The system-level variables are frequently
energy related, and they often correspond to averaged
time-domain quantities. Dynamic phasors are a uniquely useful
tool in describing mappings from system-level to micro-level
quantities (commonly denoted as “lifting” in the literature) and
back (denoted as “restriction”). In our example from power
electronics we will also use our our knowledge about the
order of switchings to properly account for ripple in microlevel
time-domain quantities to improve accuracy of the lifting
operation.
B. Structural Analysis of Energy Processing Systems
The structural analysis of energy processing systems using
tools from nonlinear dynamics (e.g., bifurcations and center
manifolds) has been used by energy engineers for two decades.
It helped elucidate important phenomena in power systems
(e.g., voltage collapse) and in power electronic systems (e.g.,
chaotic operation of PWM controlled converters). In terms of
our terminology, all applications to date have been equationbased,
and our equation-assisted framework can be seen
as complementary. Difficulties that arise when conventional,
equation-based tools are applied to practical power systems are
detailed in the introduction of [18]. The issues of proprietary
information and legacy code have been well appreciated in
power systems community even prior to the recent economic
restructuring. The downside of reliance on models in conventional
form is clearly seen in post-mortem studies of
power systems blackouts, when long model tuning is typically
needed to replicate even the coarse features of the actual
event [19], [20]. It is also difficult to distill the roots of the
network’s global behavior from such simulation studies. Power
electronic systems exhibit a very rich set of dynamic behaviors,
and their bifurcation analysis requires careful derivation
of elaborate discrete maps. The stability of a converter is
strongly dependent upon the so-called saltation matrix, the
state transition matrix relating the state just after the switching
to the one just before. Two recent reviews focusing on dc/dc
converters are given in [21], [22], and an example involving
interaction between electrical and mechanical subsystems of
an electric drive is presented in [23]. As an illustration of
bifurcation studies in power systems, we mention the voltage
stability example from [24]; the example originated in [25],
and physical justifications for values of its parameters were
discussed in [26]. A low-order example of chaotic behavior of
a single-machine infinite-bus system was presented in [27].
Our aim is to use bursts of micro-level simulations as a
magnifying lens that would detect events related to loads (e.g.,
drives stalling and disconnecting), generators (e.g., exciters
hitting limits) and protection (e.g., correlated decisions) that
are not transparent in standard models, but may prove critical
for system-level qualitative behavior.
C. Modeling and Simulation with Dynamic Phasors
The procedure that we perform to obtain dynamic phasor
models is based on the property [28] that a possibly complex
time-domain waveform x(σ) can be represented on the interval
σ ∈(t − Ts, t] using a Fourier series of the form (“lifting”):
where ωs = 2π/Ts and Xk(t) are the complex
Fourier
coefficients, which we shall also refer to as phasors. These
Fourier coefficients are functions of time since the interval
under consideration slides as a function of time. The k-th
coefficient (or k-phasor) at time t is determined by the
following averaging (“restriction”) operation:
Our approach provides dynamic models for the dominant local
Fourier series coefficients as the window of length Ts slides
over the waveforms of interest. Note that dynamic phasors can
be interpreted as outputs of linear time-invariant (LTI) filters
with causal impulse responses [29]. Given that LTI filtering
commutes with differentiation, if the waveform x(t) admits
a linear time invariant dynamic model, say
Bu(t), then its dynamic phasor Xk(t) satisfies a similar
equation, viz.,
where Uk(t) is the dynamic phasor corresponding to the
input u. Approximate dynamic phasor models can be derived
for many standard nonlinearities using describing functions.
For example, the Fourier phasor set of a product of two timedomain
variables is obtained from a discrete-time convolution
of corresponding phasor sets of each component.
Note that if the set of phasors used for reconstruction
(lifting) in (1) is small (e.g., only the dc component in our
power electronic example), then additional information (e.g.,
about the approximate magnitude and phase of the switching
ripple) may be useful in improving its accuracy. Dynamic
phasors have been used to model many dynamic processes
in energy processing - see for example [30], [6].
IV. EXAMPLE - A SWITCHED POWER ELECTRONIC CONVERTER
A. Modeling the Converter
The up/down converter, shown in Fig. 2, is a frequently used
dc/dc circuit topology with many variants. Here we consider
sampled-data (discrete-time) dynamic models that correspond
to two distinct modes of operation. In the first, the so called
continuous conduction mode, the current in the inductor L
never falls to zero (it is constrained to be non-negative by
circuit topology). This mode corresponds to larger loads (by
power), or to smaller values of load resistance R. In the
second, discontinuous conduction mode, the load resistance is
large enough so that the inductor current has zero value in the
last portion of the period. This circuit is typically controlled
by varying the relative portion of the cycle when the transistor
is on, using pulse-width modulation (PWM). Thus, in the kth
switching cycle, the inductor current increases while the
transistor is on (i.e., when the switching function q(·) = 1),
during time interval d[k]T, where T is the switching period
and d is the duty ratio (or the zero-th dynamic phasor, or
local time average of q(·)), then decreases when the transistor
is turned off (and the diode is on), and reaches zero before
the transistor is turned back on again.
For simplicity, we consider open-loop operation from control
standpoint, i.e., the duty ratio d[·] is held constant; the
addition of control action would change the final values, but
not affect the transient qualitatively (assuming stable operation
throughout). In terms of parameter values, we are using the
example from [31, Ex. 11.1, p. 255], so that the switching
frequency is 50 kHz (so T = 20μs, C = 220μF, L = 0.25mH,
vin = 12V, < vo >= −9V, and R = 2Ω. It turns out that the
boundary between the two modes of operation corresponds
to R = 25Ω. We consider a transient from discontinuous to
continuous operation when R changes 100 → 2Ω
MICRO-LEVEL MODEL OF THE UP/DOWN CONVERTER
We start by writing down the sampled-data model for the
case of continuous conduction [31] - let x = [iL vC]T so
and q[k] is the (0−1) switching function in the k-th step. To
achieve a good time resolution, T in (4) has to be at least two
orders of magnitude shorter than the duration of the switching
period, which equals Ts = 20μs.
MACRO-LEVEL MODELS OF THE UP/DOWN CONVERTER
In this simple example we have explicit models corresponding
to the two modes of operation. For the continuous
conduction, the standard macro-level, averaged model is of
the same form as (4), with the switching period Ts replacing
T and the duty ratio d[·] replacing the switching function q[·].
This model describes the dynamic evolution of the DC phasors
IL,0, VC,0; we could also add a model for the fundamental
component of the ripple IL1, VC1 that would be in the form
of (3), with A = 1/Ts(F − I), where I is the appropriate
identity matrix. This simplicity comes in part from the open-loop
nature of the model, as the addition of closed-loop control
would make for a more involved averaging procedure.
In the case of discontinuous conduction, with the averaged
capacitor voltage VC,0 selected as the only state, the sampled-data
model with sampling time Ts is nonlinear [31]
B. Equation-Free Modeling and Analysis
In this simple example it is clear that dc-averaged (“0-
th” order dynamic phasors) inductor voltage and capacitor
current are natural system-level variables. We also have the
system-level equation-based model of the form (4) that is used
to check predictions of the equation-free approach. At the
micro-level, we use the switched model (4), with a suitably
selected time-step. We focus on the load change R : 100Ω →
2Ω which corresponds to a transition from discontinuous to
continuous conduction; in Fig. 3 we show the transient -
both switched and averaged quantities are displayed. The
transient lasts approximately 250 switching periods, and the
ripple is quite pronounced for the inductor current - this will
necessitate somewhat longer reporting times τ in the equationfree
approach (and thus a larger computational effort).
In our equation-free analysis of steady-states we use two
Matlab modules - one implements the time-stepper (4) together
with the lifting and restriction operations (1)-(2), while the
other is Newton-Raphson code that includes a GMRES routine
for solving linear equations. In Fig. 4 we show the convergence
of the Newton-Raphson procedure to the (correct) steady state
values for the converter. After accounting for somewhat longer
reporting times τ that include several switching cycles (necessitated
by the switching ripple), the overall computational
savings over a direct switched simulation are of the order of
five times.
A major advantage of the equation-free framework is its
ability to incorporate approximate system-level model. Such
models (e.g., approximate Jacobians and their inverses, computed
off-line) are used to improve numerical conditioning of
problems that have to solved, for example, in the Newton-
Raphson procedure. We use (4) for this purpose, but we offset
R by 50% to make the process more realistic. In Fig. 5 we
again show the convergence of the Newton-Raphson procedure
to the steady state. The convergence is now much improved,
suggesting that computational savings of one order of magnitude
are entirely achievable in equation-assisted system-level
modeling and analysis.
As our final example, we evaluate eigenvalues (of the
explicitly
unavailable) linearization of (4) for different sampling
times, and show them in Fig. 6 together with the ones calculated
from the analytical model. We used our time-stepper in
conjunction with a matrix-free eigenvalue calculation routine
available in Matlab (command eigs). The overall agreement
in predicting these complex eigenvalues is quite good. Please
note that the RC time constant is of the order of 20 PWM
periods, so eigenvalues of linearized sampled data models for
reporting periods τ that approach the RC constant of the circuit
are of little value.
V. CONCLUSIONS
The paper outlines a procedure for system-level modeling
of
energy processing systems that is equation-free, and thus relies
on a smaller set of modeling assumptions than the conventional
system-level methods. It is in many ways complementary
to standard tools, as it has the potential to be relevant for
cases that evaded satisfactory modeling thus far. Equation-free
wrappers can be added to existing modeling tools to broaden
their applicability and to provide additional analytical features.
While the computational effort for the procedure tends to be
reasonable, it also needs to be further explored for truly large
systems.
VI. ACKNOWLEDGMENTS
This work has been supported in part by the National
Science
Foundation under Grant ECS-0601256. I thank Professor Ioannis
Kevrekidis of Princeton University for guidance regarding
the equation-free methodology
