CE 295 — Energy Systems and Control Professor Scott Moura — University of California, Berkeley

CHAPTER 2: MACHINE LEARNING

1 Overview

What is machine learning (ML)? From our perspective, it is quite simply mathematical modeling
from data. Recall Chapter 1 on mathematical modeling. Recall that we outlined seven steps as a
guide for creating mathematical models. Step 6 involves training the unknown system parameter
from experimental data. This chapter dives deeply into this specific task. Note, the mathematical
model itself may arise from either first-principles, or be a generic “black-box” format. In either case,
our objective is to “fit” the unknown parameters (a.k.a. “weights” or “coefficients”) from experimen-
tal data. In energy systems, we sometimes find ourselves rich in data yet poor in engineering
insight. These techniques provide the engineer with various tools to generate system models from
data.
Machine learning is, in fact, a convergence of several other intellectual fields that include op-
timization, statistics, control theory, and computer science. We have strategically introduced op-
timization before ML in this course, because ML can be viewed as an application of optimization.
Historically, “machine learning” was coined by IBM computer scientist Arthur Samuel in 1959.
Specifically, Arthur Samuel said the following about ML: “It gives computers the ability to learn
without being explicitly programmed.” In the 1980s – 1990s, the concept of artificial intelligence
(AI) emerged as another re-branding of ML. Interestingly, AI and ML lost popularity throughout the
2000s. The main reasons include: (i) The promises of AI and ML far exceeded the results, and
therefore academia and industry disinvested in these ideas. (ii) The lack of consistent results is
partially attributed to theoretical misunderstandings, many of which have since been resolved. (iii)
We now have the computational resources to handle very large datasets and parallelized compu-
tation. Consequently, ML and AI have re-emerged in the 2010’s as “Data Science”. Now academia
and industry are re-investing heavily in these technologies. In this chapter, you will see that ML
is essentially an application of optimization, statistics, and control theory. With this view, you will
become ML experts who understand the fundamental principles.
To this end, Chapter 2 examines the following questions:

• How do we formulate a parametric model?

• How do we identify a parametric model from data?

• Why do these algorithms work? Why do they sometimes not work?

1.1 White-box vs. Black-box Models

Inferring models from observations and studying their properties is the essential basis of science.
These models (a.k.a. hypotheses, laws of nature, theories) can be more or less formal in nature.

Revised March 30, 2018 | NOT FOR DISTRIBUTION Page 1

CE 295 — Energy Systems and Control Professor Scott Moura — University of California, Berkeley

However, they all attempt to link observations together through some pattern. Machine learning
involves the inference of a system model from input-output observations, as shown in Fig. 1. The
construction of this system often falls within three different paradigms: white-box models, black-
box models, or grey-box models.

u y
- Mathematical Model |Σ -
Input variables Output variables
Parameters |θ

Figure 1: Machine Learning involves inferring Σ and/or θ, from measurements of u, y.

• [White-box Models]: Sometimes we can formulate the mathematical model equations from
first principles (e.g. Newton’s laws, Navier-Stokes, Maxwell’s equations, etc.). This white-box
model paradigm is often preferable, since it simplifies the model identification process and
endows the model with physically relevant meaning. Moreover, these models are capable of
accurately predicting behavior outside of the training data.

• [Black-box Models]: In some applications, the underlying dynamic process is too complex
to model from first principles. Examples include human behavior, biological systems, social
systems, or large-scale systems with many intertwined dynamic phenomena. The black-
box modeling paradigm starts from measurements of the system’s inputs and outputs and
determines a mathematical relation between them without attempting to model the internals
of the system. Note that these models are generally unreliable extrapolating behavior outside
of the training data range.

• [Grey-box Models]: The aforementioned categories are end-points on a spectrum. In many

practical applications, we seek a combination of white-box and black-box models. This is the
most interesting and practical category. Grey-box models have some components that we
can model from first principles, and other components that must be fit empirically with data.

In this chapter, we are interested in tools applicable to white-box models, black-box models,
and all the shades of grey in between. Interested readers may refer to the textbook by Ljung [1]
for a thorough exposition of system identification methods for black-box models.

1.2 Online vs. Offline Algorithms

Throughout this chapter, we introduce offline and online algorithms for machine learning. Offline
(a.k.a. batch or non-recursive) algorithms operate on a fixed collection of data, asynchronous from

Revised March 30, 2018 | NOT FOR DISTRIBUTION Page 2

CE 295 — Energy Systems and Control Professor Scott Moura — University of California, Berkeley

the true physical system. These algorithms are often simpler to formulate. However, they require
large collections of data and computational power, and are not amenable to implementation in an
embedded computational system. Sometimes one wishes to “learn” the model in a continuous
manner, where it utilizes data as it arrives and then discards it. These techniques are known as
online (a.k.a. recursive) algorithms. These algorithms are typically more difficult to formulate, but
are amenable to implementation in embedded systems in synchrony with the dynamic system.

1.3 Chapter Organization

The remainder of this chapter is organized as follows:

1. Parametric Modeling

2. Gradient Algorithm

3. Least-Squares (LSQ) Algorithm

4. Nonlinear Least-Squares (NL-LSQ) Algorithm

5. Offline Techniques

6. Sensitivity Analysis

2 Regression

We begin our exposition of regression with the following example.

Example 2.1 (Linear Regression). Suppose you have collected measured data pairs (xi , yi ), for
i = 1, · · · , N where N > 6, as shown in Fig. 2. You seek to fit a fifth-order polynomial to this data,
i.e.
y = c0 + c1 x + c2 x2 + c3 x3 + c4 x4 + c5 x5 (1)

A crucial point we emphasize is that (1) is linear-in-the-parameters ci . Indeed, the output y is

nonlinear with respect to x. However, our focus is on fitting values for the ci ’s, and y is linear with
respect to each and every ci .
The goal is to determine parameters ci , i = 0, · · · , 5 that “best” fit the data in some sense. To
this end, you may compute the residual r for each data pair:

c0 + c1 x1 + c2 x21 + c3 x31 + c4 x41 + c5 x51 − y1 = r1 ,

c0 + c1 x2 + c2 x22 + c3 x32 + c4 x42 + c5 x52 − y2 = r2 ,
.. .
. = ..
c0 + c1 xN + c2 x2N + c3 x3N + c4 x4N + c5 x5N − yN = rN , (2)

Revised March 30, 2018 | NOT FOR DISTRIBUTION Page 3

CE 295 — Energy Systems and Control Professor Scott Moura — University of California, Berkeley

Figure 2: You seek to fit a fifth-order polynomial to the measured data above.

which can be arranged into matrix-vector form Ac − y = r, where

 
c0
     
1 x1 x21 x31 x41 x51 y1 r1
 

 c1 

1 x2 x22 x32 x42 x52 y2 r2
     
   c2     
 . (3)
A= , c= , y= , r=
 
.. .. .. .. .. .. .. ..

 . . . . . .

 
 c3 


 .



 .


1 xN x2N x3N x4N x5N c4  yN rN

 
c5

Now we compute an optimal fit for c in the following sense. We seek the value of c which minimizes
the squared residual

1 1 1 1 1
min krk22 = rT r = (Ac − y)T (Ac − y) = cT AT Ac − y T Ac + y T y. (4)
c 2 2 2 2 2

Note that (4) is quadratic in variable c and therefore a convex function of c. This produces a
quadratic program. In this case the problem is unconstrained. That is, our immediate objective
is to solve an unconstrained QP to fit the parameters. To this end, we set the gradient with respect
to c to zero and directly solve for the minimizer.

∂ 1
krk22 = AT Ac − AT y = 0,
∂c 2
AT Ac = AT y,

c = (AT A)−1 AT y (5)

Revised March 30, 2018 | NOT FOR DISTRIBUTION Page 4

CE 295 — Energy Systems and Control Professor Scott Moura — University of California, Berkeley

This provides a direct formula for fitting the polynomial coefficients ci , i = 0, · · · , 5 using the
measured data.
Exercise 1. Consider fitting the coefficients c1 , c2 , c3 of the following sum of radial basis functions
to data pairs (xi , yi ), i = 1, · · · , N .

2 2 2
y = c1 e−(x−0.25) + c2 e−(x−0.5) + c3 e−(x−0.75) (6)

Formulate the corresponding QP problem. Derive the explicit formula for fitting the coefficients.
Exercise 2. Repeat the same exercise for the following Fourier Series:

y = c1 sin(ωx) + c2 cos(ωx) + c3 sin(2ωx) + c4 cos(2ωx) (7)

Exercise 3 (Tikhonov or L2 regularization, a.k.a. Ridge Regression). Consider the fifth-order

polynomial regression model in (1). Suppose we seek the value of c which minimizes the squared
residual plus a so-called Tikhonov regularization term:

1 1
min krk22 + kΓck22 . (8)
c 2 2

for some matrix Γ ∈ R6×6 . Derive the QP. Solve for the minimizer of this unconstrained QP. Provide
a formula for the optimal coefficients c.

3 Maximum Likelihood Estimation

4 Parametric Modeling

Algorithms and optimization receive all the glory in model identification. However, skillful paramet-
ric modeling gets grossly under-appreciated. It is a necessary initial step for any model identifica-
tion study. A cleverly formulated parametric model can render the algorithms nearly trivial. This
section discusses various canonical formulations for parametric models.
In Chapter 1, we discussed the derivation of models from first principles. This typically renders
ODEs that can be formulated into state space form as a series of first-order ODEs

ẋ(t) = f (x(t), u(t); θ) (9)

where x(t) ∈ Rn×1 , u(t) ∈ Rp×1 , t represents continuous time, and θ ∈ Rp×1 are a set of constant
parameters to be identified. In practice, data is sampled in discrete-time instances. Consequently
we shall also consider discrete-time models in this chapter, which take the form

x(k + 1) = f (x(k), u(k); θ) (10)

Revised March 30, 2018 | NOT FOR DISTRIBUTION Page 5

CE 295 — Energy Systems and Control Professor Scott Moura — University of California, Berkeley

In this so-called “difference equation model”, k is related to t according to t = k∆t, where ∆t is

the sampling time.

4.1 Linear Parametric Models

Consider the discrete-time analog of an n−th order differential equation given by

y(k + 1) + a0 y(k) + a1 y(k − 1) + . . . + an y(k − n) = b0 u(k) + b1 u(k − 1) + . . . + bn u(k − n) (11)

where ai , bi , for i = 0, 1, . . . , n are scalar parameters. In the stochastic process literature [1], this
is known as the famous ARX model (auto-regressive with exogenous input u()), and is a common
black-box model. If we lump all the ai , bi parameters into the vector

θ = [b0 , b1 , . . . , bn , a0 , a1 , . . . , an ]T (12)

and all the input/output (I/O) signals into the vector

φ(k) = [u(k), u(k − 1), . . . , u(k − n), −y(k), −y(k − 1), . . . , −y(k − n)]T (13)

then we may re-write (11) into matrix-vector form as

y(k + 1) = θT φ(k). (14)

Equation (14) is linear in θ which, as we show in subsequent sections, is crucial for designing
parameter identifiers to estimate θ from the measurements of y(k + 1) and φ(k).

4.2 Nonlinear Parametric Models

It is always beneficial to formulate the model into the form (14) when possible. However, some-
times the model is not linear in the parameters. In these cases, the set of available algorithms
decreases dramatically and their properties generally weaken. However, these nonlinear cases
do arise in practice. Generally, a nonlinear in the parameters model takes the form of

y(k + 1) = f (φ(k); θ) (15)

where φ(k) is an appropriately defined vector of I/O data, y(k + 1) is the output data at time-step
k + 1, θ is the unknown parameter vector, and f (·, ·) is a nonlinear function that maps I/O data and
parameters to a scalar measured output. Artificial neural networks are an example of a nonlinear-
in-the-parameters model. In this case, we must utilize iterative optimization methods to estimate
the value of parameter vector θ that best matches the data and model.

Revised March 30, 2018 | NOT FOR DISTRIBUTION Page 6

CE 295 — Energy Systems and Control Professor Scott Moura — University of California, Berkeley

x1
@
@
@ w1
@ #
x2 w2 @R
@
- f (·)
y -
"!
w3

x3

Figure 4: Simple single neuron network with acti-

Figure 3: An artificial neural network (ANN) is a
vation function f (·).
weighted and directed graph, where the nodes are
called “neurons”.

Example 4.1 (Artificial Neural Networks (ANN)). Artificial neural networks (ANNs) were first cre-
ated by Warren McCulloch and Walter Pitts in 1943 as a computational model for the human brain.
Recently, ANNs have re-emerged within artificial intelligence due to (i) the availability of enhanced
computational power, and (ii) resolutions to theoretical stability issues.
A neural network is composed of a weighted directed graph, as shown in Fig. 3. Each node of
the graph (a.k.a. neuron) is “activated” by inputs. The node will perform some simple computations
on these inputs, using a so-called “activation function”, and then pass its result along the network.
More specifically, each node computes a weighted sum of its inputs and then processes this sum
with its activation function, i.e. f ( i wi xi ) or written in vector notation as f (wT x).
P

Consider the trivially simple single neuron network shown in Fig. 4. Denote the input weights
by wi . There are several canonical activation functions.
Linear-in-the-parameters Model: Consider the affine activation function

f (wT x) = a · wT x + b, (16)

where a, b ∈ R. This renders the parametric input-output model: y = a · wT x + b. Now, suppose

we define the parameter vector to be:

θ = [a, b, w1 , w2 , · · · , wn ]T . (17)

where θ ∈ Rn+2 . This is intuitive. We just collect each scalar parameter into a vector θ. However,
this yields a parametric input-output model which is nonlinear in θ. To see this, let θi be the ith

Revised March 30, 2018 | NOT FOR DISTRIBUTION Page 7

CE 295 — Energy Systems and Control Professor Scott Moura — University of California, Berkeley

component of θ. Then we have

n−1
X
y = θ1 · θj+2 xj + θ2 (18)
j=1

where θ1 · θj+2 is a nonlinear term.

An alternative choice (and better w.r.t. learning convergence) would be to define the parameter
vector to be:
ϑ = [a · w1 , a · w2 , · · · , a · wn , b]T . (19)

where ϑ ∈ Rn+1 . Let ϑi be the ith component of ϑ. Then we have

n
X
y= ϑi xi + ϑn+1 (20)
i=1

This parametric model is linear w.r.t. ϑ. Moreover, it has less parameters n + 1 versus n + 2 for θ.
Nonlinear-in-the-parameters Model: Two other common examples are the sigmoid and tanh
activation functions, given respectively by

1
f1 (wT x) = , f2 (wT x) = tanh(wT x). (21)
1 + exp(wT x)

These two activation functions produce parametric input-output models given by:

1
y= , y = tanh(wT x), (22)
1 + exp(wT x)

which are both nonlinear in parameter vector w.

Example 4.2 (Flywheel Energy Storage). Consider the flywheel energy storage example from
CH1. Recall that the equation of motion is governed by Euler’s rotation equation:

I ω̇(t) = −bω(t) + T (t), (23)

˙ denotes the derivative with respect to time. Note that ω(t) is the state variable, T (t) is a
where (·)
controllable input torque, and I, b are physical parameters representing the rotational inertia and
friction coefficient, respectively. Suppose we can measure ω(t), T (t). Our ultimate goal is to learn
the unknown parameters. To this end, we must first derive a parametric model. In this example,
we will derive both linear and nonlinear-in-the-parameters models.
d ω(t+∆t)−ω(t)
Suppose we approximate the time-derivative using forward Euler’s method: dt ω(t) ≈ ∆t .
Furthermore, we define discrete-time index k such that t = k · ∆t. Then we can re-write (23) as

Revised March 30, 2018 | NOT FOR DISTRIBUTION Page 8

CE 295 — Energy Systems and Control Professor Scott Moura — University of California, Berkeley

discrete time difference equation

ω(k + 1) − ω(k)
I· = −bω(k) + T (k), (24)
∆t  
b∆t ∆t
ω(k + 1) = 1 − ω(k) + T (k) (25)
I I

Linear-in-the-parameters Model: The flywheel dynamics are now in the form of ω(k + 1) = θT φ in
(14), where the parameter vector θ and regressor signal φ are:
 T
b∆t ∆t
θ = 1− , , φ = [ω(k), T (k)]T . (26)
I I

Interesting, if the sampling time ∆t is known (a very mild assumption), then we can back-calculate
b∆t ∆t ∆t 1−θ1
I, b from θ1 = 1 − I , θ2 = I by the following formulae: I = θ2 , b = θ2 . Note, it is NOT
always possible to back-calculate the physical parameters from the linear parameter vector θ. A
necessary (but not sufficient) condition for back-calculating physical parameters from θ is that the
number of physical parameters equals the dimension of θ.
Nonlinear-in-the-parameters Model: Alternatively, we might initially define the parameter vector to
correspond to the physical parameters as follows: θ = [I, b]T . In this case, we have a nonlinear-
in-the-parameters model given by:

ω(k + 1) = f (ω(k), T (k); θ) (27)

Loosely speaking, it is MUCH more difficult to identify parameters for the nonlinear-in-the-parameters
model above. Technically speaking, there are no general algorithms that are provably convergent.
Hence, linear-in-the-parameters models are always preferred.

4.3 Identifiability

A logical question to ask is: “Is it always possible to determine the model parameters by observing
a sufficient amount of data? If not, then is there a test to determine when it is possible to determine
the model parameters?” In model identification this existence question is called identifiability. We
now present several notions of identifiability:

• [Identifiable]: A model is said to be identifiable if it is possible to uniquely learn the param-

eters by observing a sufficient amount of data.

• [Partially Identifiable]: A model is said to be partially identifiable if it is possible to uniquely

learn a sub-set of the parameters by observing a sufficient amount of data.

Note that identifiability and partial identifiability depend on your data. When people say “I have
excellent data quality,” they are referring (perhaps ignorantly) to this concept. Given these con-

Revised March 30, 2018 | NOT FOR DISTRIBUTION Page 9

CE 295 — Energy Systems and Control Professor Scott Moura — University of California, Berkeley

cepts, a simple test to determine identifiability would be extremely convenient. We shall discuss
two tests in this chapter. The first test given in Section 5.3 is specific to models that are linear-in-
the-parameters, and is known as the persistence of excitation condition. The second test given in
Section 9 is applicable to nonlinear-in-the-parameters models, and is called sensitivity analysis.

5 Gradient Algorithm

5.1 Scalar Example

Consider the scalar input-output system given by the algebraic equation

y(t) = θu(t) (28)

where y(t) is a scalar output, u(t) is a scalar input, and θ is an unknown scalar parameter. As-
suming that u(t), y(t) are measured, we seek an estimate of θ at each time t. We denote the
estimated version of the parameter vector as θ̂. If the measurements of y, u were noise free, one
could simply calculate θ̂(t) as
y(t)
θ̂ = (29)
u(t)
whenever u(t) 6= 0. The division (29), however, may not be desirable because u(t) may assume
values arbitrarily close to zero. Furthermore, the effect of noise on the measurement of u, y may
lead to an erroneous estimate of θ. The noise and divide-by-zero effects in (29) may be reduced
by using various other nonrecursive or offline methods, especially when θ is a constant for all t.
In our case, we are interested in a recursive or on-line method to generate θ̂(t). We are looking
for a differential equation, which (i) depends on signals that are measured, (ii) whose solution is
θ̂(t), (iii) has equilibrium point θ̂eq = θ, and (iv) this equilibrium point is asymptotically stable. The
procedure for developing such a differential equation is given below.
Using θ̂(t) as the estimate of θ at time t, we generate the estimated or predicted value ŷ(t) of
the output y(t) as
ŷ(t) = θ̂(t)u(t) (30)

The prediction or estimation error (t) is given by

(t) = y(t) − ŷ(t) = y(t) − θ̂(t)u(t) (31)

Now define the parameter estimation error as θ̃(t) = θ − θ̂(t), where the error (t) can be re-written
as
(t) = θu(t) − ŷ(t) = θ̃(t)u(t) (32)

The differential equation for generating θ̂(t) is developed by minimizing various cost criteria of 

Revised March 30, 2018 | NOT FOR DISTRIBUTION Page 10

CE 295 — Energy Systems and Control Professor Scott Moura — University of California, Berkeley

with respect to θ̂. For this example, we use the simple cost criterion

1 1 2
J(θ̂) = 2 = y − θ̂u (33)
2 2

It is easy to see that, for each time t, the cost function J(θ̂) is convex over all R; therefore, any local
d
minimum of J is also global and satisfies J(θ̂) = 0. Consequently, we may utilize a gradient-
dθ̂
scheme (a.k.a. steepest descent) to generate an update law for θ̂(t), as follows

˙ d
θ̂(t) = −γ J(θ̂) = γ(y − θ̂u)u = γu, θ̂(0) = θ̂0 (34)
dθ̂

where γ > 0 is a scaling constant. In the literature, the differential equation (34) is referred to as
the update law or identification algorithm for updating θ̂(t) to estimate θ online.
Stability analysis helps us show that the solution of (34) indeed converges to the true value of
θ, under appropriate conditions. Stability is analyzed by re-writing (34) in terms of the parameter
error θ̃ = θ − θ̂, i.e.
˙ ˙
θ̃ = θ̇ − θ̂ = θ̇ − γu (35)

Because θ is constant, i.e. θ̇ = 0, and  = θ̃u from (32), we have

˙
θ̃(t) = −γu2 (t)θ̃(t) (36)

Note that the coefficient −γu2 (t) is always non-positive, which from Chapter 1 we learned proves
that the parameter error dynamics are marginally stable or stable in the sense of Lyapunov. If we
add the sufficient condition that u(t) 6= 0 uniformly in t, then we obtain the stronger condition that
the estimation error is asymptotically stable. In other words, θ̂(t) → θ as t → ∞.

5.2 Vector Case

Next we present the vector case of the gradient algorithm. Consider a linear parametric model of
the form
z(t) = θT φ(t) (37)

where z(t) is a scalar measurement, θ is a vector of unknown parameters, and φ(t) is a vector
of I/O data. Note this is the canonical linear in the parameters model from Section 4.1. The
corresponding prediction or estimation error (t) is given by

(t) = z(t) − ẑ(t) = z(t) − θ̂T (t)φ(t) (38)

Revised March 30, 2018 | NOT FOR DISTRIBUTION Page 11

CE 295 — Energy Systems and Control Professor Scott Moura — University of California, Berkeley

C
R1
_ V +
c

+
OCV(z) R2 V
_
I

Figure 5: OCV-R-RC Equivalent circuit model of a battery. Includes an open circuit voltage in series with a
resistance, in series with a resistor-capacitor pair.

Now define the parameter estimation error as θ̃(t) = θ − θ̂(t), where the error (t) can be re-written
as
(t) = θT φ(t) − ẑ(t) = θ̃T (t)φ(t) (39)

The differential equation for generating θ̂(t) is developed by minimizing various cost criteria of 
with respect to θ̂. For this example, we use the simple cost criterion

1 1 2
J(θ̂) = 2 = z − θ̂T φ (40)
2 2

It is easy to see that, for each time t, the cost function J(θ̂) is convex over all R; therefore, any local
d
minimum of J is also global and satisfies J(θ̂) = 0. Consequently, we may utilize a gradient-
dθ̂
scheme (aka steepest descent) to generate an update law for θ̂(t), as follows

˙ d dJ d
θ̂(t) = −Γ J(θ̂) = −Γ = Γφ(t)(t), θ̂(0) = θ̂0 (41)
dθ̂ d dθ̂

where Γ = diag(γ1 , γ2 , · · · , γn )  0 is a diagonal, positive-definite matrix.

Example 5.1 (Battery Voltage). In HW1, you studied an equivalent circuit model shown in Fig. 5,
where the battery voltage output function is,

V (t) = OCV (z(t)) + Vc (t) + R1 I(t) (42)

The state variables are state-of-charge z(t) and capacitor voltage Vc (t). The input is current I(t).
In this equation, R is a parameter and function OCV (z) can be parameterized with a cubic poly-
nomial: OCV (z) = p0 + p1 z + p2 z 2 + p3 z 3 . Substituting,

V (t) = p0 + p1 z(t) + p2 z 2 (t) + p3 z 3 (t) + Vc (t) + R1 I(t) (43)

Our goal is to use the gradient algorithm to recursively learn parameters p0 , p1 , p2 , p3 , R1 in the
model (43), given streaming data V (t), z(t), Vc (t), I(t).

Revised March 30, 2018 | NOT FOR DISTRIBUTION Page 12

CE 295 — Energy Systems and Control Professor Scott Moura — University of California, Berkeley

First, we must define a parametric model. A linear-in-the-parameters model z(t) = θT φ(t) can
be defined as follows:
T
θ = [p0 , p1 , p2 , p3 , R1 ]T , φ(t) = 1, z(t), z 2 (t), z 3 (t), I(t)

z(t) = V (t) − Vc (t), (44)

Second, we apply the gradient algorithm given by

˙
θ̂(t) = Γφ(t)(t), θ̂(0) = θ̂0 (45)

where Γ = diag(γ1 , γ2 , γ3 , γ4 , γ5 )  0 is a diagonal, positive-definite matrix. Think of Γ as user-

selected parameter that dictates the “learning rate”. Also,  = z(t) − θ̂T (t)φ(t) is the estimation
error and θ̂0 is an initial guess for uncertain parameters θ.
Third, we implement the gradient algorithm with parameters Γ =diag(10−1 , 101 , 101 , 101 , 10−1 ),
θ̂0 = [3.5, 1, 0, 0, 0]T . The results are provided in Fig. 6, where the exit estimates are θ̂tf inal =
[3.504, 1.0185, −0.044, −0.037, 0.050]T . Recall the true parameter values are θ = [3.4707, 1.6112,
−2.6287, 1.7175, 0.05]T . Observe the predicted voltage V̂ (t) = p̂0 (t) + p̂1 (t)z(t) + p̂2 (t)z 2 (t) +
p̂3 (t)z 3 (t)+Vc (t)+ R̂1 (t)I(t) converges to the true voltage V (t). However, the parameter estimates
θ̂(t) do not all converge to the true values θ. This may or may not be an issue, depending on
your particular objective. Is it critical to estimate the TRUE physical parameters, or simply obtain
accurate predictions of battery voltage? In this example, the model parameters are only partially
identifiable from this model and data. We discuss identifiability next.

5.3 Identifiability Test

Similar to the scalar case, we require an extra technical condition to ensure parameter vector θ̂(t)
converges asymptotically to the true parameter θ, given sufficient input-output data. As before, we
can write the parameter estimation error dynamics as

˙ ˙
θ̃(t) = θ̇ − θ̂ = −Γφ = −Γφ(θT φ − θ̂T φ) = −Γφ(φT θ − φT θ̂),
= −ΓφφT θ̃,
= −ΓA(t)θ̃(t) (46)

where Γ is a diagonal positive-definite matrix and A(t) = φ(t)φT (t). Recalling our notion of asymp-
totic stability for an autonomous linear time-invariant (LTI) systems of the form ẋ(t) = ALT I x(t)
from Chapter 1, we require that ALT I have eigenvalues that all have strictly negative real parts. In
the present situation, we have a linear time-varying system, which requires a slightly more techni-
cal condition, known as Persistent Excitation. Intuitively, it encapsulates the fact that the data must
be sufficiently rich to excite the dynamics such that the true parameters can be inferred.

Definition 5.1. (Persistent Excitation (PE)) The (bounded) regressor signal φ(t) is called PE with

Revised March 30, 2018 | NOT FOR DISTRIBUTION Page 13

 CE 295 — Energy Systems and Control Professor Scott Moura — University of California, Berkeley

0 4.0
V(t)
3.9
1 V̂(t)
3.8
2
I(t) [amps]

3.7

Voltage
3 3.6
3.5
4
3.4
50 100 200 300 400 500 600 3.30 100 200 300 400 500 600
4.0 0.10
3.5 R̂1 (t)
3.0 0.08
2.5 p̂0 (t)
Parameters

Parameters
2.0 p̂1 (t) 0.06
1.5 p̂2 (t)
0.04
1.0 p̂3 (t)
0.5 0.02
0.0
0.50 100 200 300 400 500 600 0.000 100 200 300 400 500 600
Time [sec] Time [sec]

Figure 6: Parameter estimates of HW1 voltage model, using the gradient algorithm.

a level of excitation α if there exists some constant T such that

Z t+T
1
φ(τ )φT (τ )dτ ≥ αI, ∀t ≥ 0 (47)
T t

If the regressor signal φ(t) is PE, then θ̂(t) → θ as t → ∞. Note that a sufficient condition to satisfy
the PE definition (47), which can be checked numerically, is
Z t+T 
T
λmin φ(τ )φ (τ )dτ >0 (48)
t

where λmin (·) signifies the minimum eigenvalue of (·).

Remark 5.1. Note that the PE condition means A(t) in the dynamic system (46) must be “positive”
in some sense. Specifically, the integrated value of A(t) over some time interval of size T must be
positive definite. The “more” positive-definite the left-hand-side of (47) is, the faster the estimate
θ̂(t) converges to the true value θ. This convergence speed can be characterized by the level of
excitation, α. We summarize this important result in the following theorem.

Theorem 5.1 (Gradient Algorithm Parameter Convergence, adopted from Corollary 4.3.1 of [2]).
˙
Consider parametric model z(t) = θ̂T (t)φ(t) and update law θ̂(t) = Γφ(t)(t) with Γ diagonal and
positive definite, as defined in (37)-(41). If φ(t) is PE according to Definition 5.1 and φ(t), φ̇(t) are

Revised March 30, 2018 | NOT FOR DISTRIBUTION Page 14

CE 295 — Energy Systems and Control Professor Scott Moura — University of California, Berkeley

bounded (i.e. stay away from infinity), then θ̂(t) → θ exponentially fast.

Proof. The proof is long and pedantic for an applied course on Energy Systems and Control. We
refer interested readers to Section 4.8 of [2].

The practical utility of this theorem is that checking if φ(t) is PE is the key step to determining
if the gradient algorithm converges exponentially fast to true parameter vector θ.

Remark 5.2. For non-identifiable models and data, the PE level is - in practice - never exactly zero.
In fact, it is usually some very small number close to zero, e.g. 10−6 . This leads to the natural
question: Does PE practically enable us to determine identifiability or non-identifiability? As an
absolute numeric value, the PE level is not very useful. However, the PE level metric allows us to
compare different data sets or different models to use for model identification.

Remark 5.3. As energy systems engineers, we can often design our input-output data. The
consequence of this is that data can be selected such that φ is PE. Checking this PE property
a priori provides a certificate that the data will indeed render estimates that converge to the true
values. More interestingly, we can seek to optimize φ(t) such that the left-hand-side of (48) is
maximum, thereby producing the richest possible data for model identification.

Example 5.2 (Battery Voltage). Consider again the equivalent circuit model shown from Ex. 5.1.
The parametric model is given by z = θT φ, where
T
θ = [p0 , p1 , p2 , p3 , R1 ]T , φ(t) = 1, z(t), z 2 (t), z 3 (t), I(t)

z(t) = V (t) − Vc (t), (49)

Using the -5A pulsed current input in Fig. 6, the PE level can be calculated via (48) to be: 10−6 .
For comparison, let us consider a linear OCV function, which gives the parametric model z = ϑT ϕ:

z(t) = V (t) − Vc (t), ϑ = [p0 , p1 , R1 ]T , ϕ(t) = [1, z(t), I(t)]T (50)

This parametric model has only 3 parameters, whereas the previous had 5 parameters. Com-
puting the PE level via (48) yields: 0.013. Consequently, the smaller parameterization has strong
identifiability properties, despite using a linear OCV function instead of the true cubic function.
We might also consider different data. Suppose we apply a constant -1A current input to the
battery model. The PE level for both parametric models yields a trivially small result: 10−16 . This
implies the pulsed current is better for model identification than the constant current. The results
are summarized in Table 1.

6 Recursive Least Squares Algorithm

Least-squares is an old method dating back to Gauss in the eighteenth century where he used it
to determine the orbit of planets. The basic idea behind the least-squares is fitting a mathematical

Revised March 30, 2018 | NOT FOR DISTRIBUTION Page 15

CE 295 — Energy Systems and Control Professor Scott Moura — University of California, Berkeley

Table 1: PE Level for Equivalent Circuit Battery Model

Data -5A pulse -1A constant

θ 10−6 10−16
ϑ 0.013 10−16

model to a sequence of observed data by minimizing the sum of the squares of the difference
between the observed and computed data. In doing so, any noise or inaccuracies in the observed
data are expected to have less effect on the accuracy of the mathematical model. Next we derive a
recursive version of the famous least squares algorithm, known as recursive least squares (RLS).
The method is simple to apply and analyze for linear-in-the-parameters models of the form

z(t) = θT φ(t) (51)

As before, consider the estimate ẑ of z and the estimation error are generated as

ẑ(t) = θ̂T (t)φ(t), (t) = z(t) − θ̂T (t)φ(t) (52)

In this case, we are interested in minimizing a cost function which integrates the squared error
over time, as follows
Z t Z th
1 2 1 i2
J(θ̂) = (τ ) dτ = z(τ ) − θ̂T (t)φ(τ ) dτ (53)
2 0 2 0

The cost J(θ̂) penalizes all the past errors from τ = 0 to t that are due to θ̂(t) 6= θ. Note the
conceptual difference between this cost function, and the cost function (40) used for the gradient
algorithm. Specifically, the gradient algorithm minimizes an instantaneous cost, whereas RLS
minimizes an integrated cost. This is a crucial point, conceptually. Since J(θ̂) is a convex function
over R at each time t, its minimizer θ̂ satisfies
Z t Z th
d dJ d T
i
0 = J(θ̂) = =− (τ )φ (τ )dτ = θ̂T (t)φ(τ ) − z(τ ) φT (τ )dτ, (54)
dθ̂ d dθ̂ 0 0
Z t Z t
0 = θ̂T (t) φ(τ )φT (τ )dτ − z(τ )φT (τ )dτ (55)
0 0

for any given time t, where dJ(θ̂)/dθ̂ ∈ R1×n is defined as a row vector. Solving for θ̂(t) results in
Z t
θ̂(t) = P (t) z(τ )φ(τ )dτ (56)
0

Revised March 30, 2018 | NOT FOR DISTRIBUTION Page 16

CE 295 — Energy Systems and Control Professor Scott Moura — University of California, Berkeley

where −1
Z t
T
P (t) = φ(τ )φ (τ )dτ (57)
0

provided the inverse exists. This is the celebrated offline (a.k.a. batch or non-recursive) least
squares algorithm.
We can derive a recursive method using the identity

d d
P P −1 = Ṗ P −1 + P P −1 = 0 (58)
dt dt

to show that P (t) satisfies the differential equation

Ṗ (t) = −P (t)φ(t)φT (t)P (t), P (0) = P0 (59)

where P0 = diag(p01 , p02 , · · · , p0n )  0 is diagonal and positive definite. The implication of (59) is
the calculation of the inverse in (57) can be avoided by generating P as the solution to the matrix
differential equation (59). Similarly, differentiating θ̂(t) w.r.t. t in (56) and using (59) we can derive
the relation
˙
θ̂(t) = P (t)φ(t)(t) (60)

To summarize, this produces the continuous-time recursive or online least-squares algorithm,

˙
θ̂(t) = P (t)φ(t)(t), θ̂(0) = θ̂0 , (61)
T
Ṗ (t) = −P (t)φ(t)φ (t)P (t), P (0) = P0 (62)

Exercise 4. Derive (59) by substituting (57) into (58) and solving for Ṗ . Next, derive (60) by
differentiating θ̂(t) w.r.t. t in (56) and using (59).

For recursive least squares (RLS), one can also prove convergence to the true parameter
values. We summarize this result with the following theorem.

Theorem 6.1 (RLS Parameter Convergence, adopted from Corollary 4.3.4 of [2]). Consider para-
˙
metric model z(t) = θ̂T (t)φ(t) and RLS update law θ̂(t) = P (t)φ(t)(t), Ṗ (t) = −P (t)φ(t)φT (t)P (t)
with P (0) = P0 diagonal and positive definite, as defined in (37) and (61)-(62). If φ(t) is PE ac-
cording to Definition 5.1 and φ(t) is bounded (i.e. stay away from infinity), then θ̂(t) → θ as t → ∞.

Proof. The proof is long and pedantic for an applied course on Energy Systems and Control. We
refer interested readers to Section 4.8 of [2].

The practical utility of this theorem is that checking if φ(t) is PE is the key step to determining
if RLS converges asymptotically fast to true parameter vector θ.

Revised March 30, 2018 | NOT FOR DISTRIBUTION Page 17

CE 295 — Energy Systems and Control Professor Scott Moura — University of California, Berkeley

Example 6.1 (Battery Voltage). Consider once again the equivalent circuit model shown from Ex.
5.1. The parametric model is given by z = θT φ, where
T
θ = [p0 , p1 , p2 , p3 , R1 ]T , φ(t) = 1, z(t), z 2 (t), z 3 (t), I(t)

z(t) = V (t) − Vc (t), (63)

We apply the RLS algorithm given by

˙
θ̂(t) = P (t)φ(t)(t), θ̂(0) = θ̂0 , (64)
Ṗ (t) = −P (t)φ(t)φT (t)P (t), P (0) = P0 (65)

where P0 = diag(p01 , p02 , p03 , p04 , p05 )  0 is a diagonal, positive-definite matrix. Think of P0 as
the initial “learning rate”, which gets automatically scaled by 65. Also,  = z(t) − θ̂T (t)φ(t) is the
estimation error and θ̂0 is an initial guess for uncertain parameters θ.
Next, we implement the RLS algorithm with parameters P0 =diag(10−1 , 101 , 101 , 101 , 10−1 ),
θ̂0 = [3.5, 1, 0, 0, 0]T . The results are provided in Fig. 7, where the exit estimates are θ̂tf inal =
[3.525, 0.898, −0.370, −0.246, 0.048]T . Recall the true parameter values are θ = [3.4707, 1.6112,
−2.6287, 1.7175, 0.05]T . Also, note that some parameters estimated by RLS are slightly different
than the exit estimates from the gradient law in Ex 5.1. Observe the predicted voltage V̂ (t) =
p̂0 (t) + p̂1 (t)z(t) + p̂2 (t)z 2 (t) + p̂3 (t)z 3 (t) + Vc (t) + R̂1 (t)I(t) converges to the true voltage V (t),
but at a slower rate than the gradient algorithm. This effect is intuitive, since RLS penalizes
the integrated square error as opposed to the instantaneous square error. In other words, the
evolution of the parameter updates θ̂(t) is less sensitive to prediction errors (t) with RLS than
with the gradient algorithm.

Algorithm Embellishments and Variations

The version of least-squares introduced above is known as “pure” least squares, to differentiate it
from the various embellished versions. For a comprehensive coverage of these various adaptive
laws, please consult [2]. These alternate versions impose constraints on the estimated parameters
(projection), weight recent data more heavily than older data (forgetting factors), and reset the
covariance (covariance resetting). A summary of various adaptive laws is given in the Appendix,
for interested readers.

Exercise 5. A popular embellishment of the least squares cost function in (53) is to add a penalty
on the parameter magnitude. These penalties ameliorate ill-conditioned problems, and are gen-
erally called “regularization” in mathematics and statistics. For example, consider the integrated
square error with weighted L2 penalty:
Z t Z th
1 1 1 i2 1
J(θ̂) =  (τ )dτ + θ̂T Γθ̂ =
2
z(τ ) − θ̂T (t)φ(τ ) dτ + θ̂T (t)Γθ̂(t) (66)
2 0 2 2 0 2

Revised March 30, 2018 | NOT FOR DISTRIBUTION Page 18

 CE 295 — Energy Systems and Control Professor Scott Moura — University of California, Berkeley

0 4.0
V(t)
3.9
1 V̂(t)
3.8
2
I(t) [amps]

3.7

Voltage
3 3.6
3.5
4
3.4
50 100 200 300 400 500 600 3.30 100 200 300 400 500 600
4.0 0.10
3.5 R̂ 1 (t)
3.0 0.08
2.5 p̂ 0 (t)
Parameters

Parameters
2.0 p̂ 1 (t) 0.06
1.5 p̂ 2 (t)
0.04
1.0 p̂ 3 (t)
0.5 0.02
0.0
0.50 100 200 300 400 500 600 0.000 100 200 300 400 500 600
Time [sec] Time [sec]

Figure 7: Parameter estimates of HW1 voltage model, using recursive least squares.

where Γ is a diagonal positive definite weighting matrix. This cost function with L2 penalty is known
as Tikhonov regularization or “ridge regression” in the statistics literature. For this exercise, derive
the updates laws for θ̂(t) and P (t) for the cost function in (66). Your answer will be analogous to
(61)-(62). HINT: The answer is exactly (61)-(62), except the initial condition for P (t) is P (0) = Γ−1 .

7 Nonlinear Least Squares Algorithm

In practice, we often encounter nonlinear-in-the-parameters models, such as (15). Examples of

nonlinear-in-the-parameters models include activation energy E in the Arrhenius law k = Ae−E/RT ,
2
shape parameter ε in the radial basis function φ(r) = e−(εr) , and frequency ω in the harmonic
function sin(ωt) + cos(ωt). In fact, most neural network models (including the recently popular
class of “deep learning networks”) are characterized by nonlinear-in-the-parameter models. Non-
linear least squares is by far the most popular algorithm for handling these models. The key steps
of this algorithm are:

1. Re-express the nonlinear model in terms of the parameter estimation error θ̃ = θ − θ̂.

2. Perform a Taylor series expansion around θ̃ = 0.

3. Define a linearized parametric model based upon a first-order Taylor series approximation

Revised March 30, 2018 | NOT FOR DISTRIBUTION Page 19

CE 295 — Energy Systems and Control Professor Scott Moura — University of California, Berkeley

Figure 8: Fuel cell polarization curve from [4].

4. Apply least squares, using a nonlinear estimation error.

We demonstrate this algorithm by example.

Example 7.1 (Fuel Cell Polarization Curve). Figure 8 presents the output voltage of a fuel cell, as
a function of the current density [A/cm2 ]. This is known as the fuel cell polarization curve. The fuel
cell voltage is calculated using a combination of physical and empirical relationships, and is given
by [3]
Vf c = E − Vact − Vohm − Vconc (67)

where E is the open circuit voltage, and Vact , Vohm , Vconc are activation, ohmic, and concen-
tration over voltages, which represent losses due to various physical or chemical factors. For
pedagogical purposes, we assume Vconc is negligible. The remaining loss terms are given by

Vact = va 1 − e−c·i ,


(68)
Vohm = Ri (69)

where E, va , c, R are unknown parameters to be identified. The problematic parameter is c, since

this parameter is nonlinear in the expression. Step 1: We redefine the parameters as θ1 = E, θ2 =
va , θ3 = c, θ4 = R and re-write the polarization curve equation in terms of the parameter error as
h i
Vf c = (θ̃1 + θ̂1 ) − (θ̃2 + θ̂2 ) 1 − e−(θ̃3 +θ̂3 )i − (θ̃4 + θ̂4 )i (70)

Revised March 30, 2018 | NOT FOR DISTRIBUTION Page 20

CE 295 — Energy Systems and Control Professor Scott Moura — University of California, Berkeley

Step 2: Next we perform a Taylor series expansion around θ̃ = 0,

h i
Vf c = θ̂1 − θ̂2 1 − e−θ̂3 i − θ̂4 i
h i
+θ̃1 − 1 − e−θ̂3 i · θ̃2 − θ̂2 ie−θ̂3 i · θ̃3 − i · θ̃4 + (higher order terms) (71)

Step 3: Truncating and re-arranging terms, we arrive at the linearized parametric model

enl = θ̃T φ (72)

where
h i
enl = Vf c − θ̂1 + θ̂2 1 − e−θ̂3 i + θ̂4 i, (73)
 
1
 −θ̂3 i 
 e − 1 
φ =  −θ̂ ie−θ̂3 i 
 (74)
 2 
−i

and both enl and φ are functions of the current parameter estimate θ̂. Step 4: is to apply the least
squares parameter update law as

˙
θ̂ = P enl φ, θ̂(0) = θ̂0 (75)
T
Ṗ = −P φφ P, P (0) = P0 = P0T >0 (76)

This algorithm has local stability properties. That is, when the parameter estimate θ̂ is near its true
value, then the linear approximation is reasonably accurate and θ̂ converges asymptotically to the
true value θ under an appropriate persistence of excitation condition for the regressor φ.

Exercise 6 (Photovoltaic Array Learning). Consider a photovoltaic (PV) array. Figure 9 [top]
presents the characteristic nonlinear relationship between the array’s voltage and current. Mul-
tiplying the current by voltage yields the concave relationship between array current and output
power in Fig. 9 [bottom]. Maximum power point tracking (MPPT) algorithms seek to maintain
the PV at maximum output power, in the face of varying ambient temperatures and incident solar
irradiation (see [5] and references therein).
In this exercise, you seek to learn the characteristic nonlinear relationship in Figure 9, using
current I(t) and voltage V (t) data. Since these curves shift with varying temperature and irra-
diation levels, you seek to learn this relationship in real-time. Derive a nonlinear least squares
algorithm to identify θ = [θ1 , θ2 , θ3 , θ4 , θ5 ]T using the following nonlinear-in-the-parameters model:
h i
θ1 − θ2 e(θ3 V (t)+θ4 I(t)) − 1 − θ5 [θ3 V (t) + θ4 I(t)] = 0 (77)

Revised March 30, 2018 | NOT FOR DISTRIBUTION Page 21

CE 295 — Energy Systems and Control Professor Scott Moura — University of California, Berkeley
(a) (b)
3 200 W/m2 3 0°C
2
400 W/m 20°C
2
2.5 (a) 600 W/m 2.5 (b) °
2 40
0°CC
3 200
800 W/m
W/m2 3 °
2 400
2 2 60
20°C
1000W/m C
Current [A]

Current [A]
W/m2
2
2.5 600 W/m 2.5
1.5 1.5 40°C
800 W/m2
2 2 60°C
Current [A]

Current [A]
1 1000 W/m2 1
1.5 1.5
0.5 0.5
1 1
0 0
0 5 10 15 20 25 0 5 10 15 20 25
0.5 0.5

0 0
0 5 10 15 20 25 0 5 10 15 20 25
50 50

40 40
50 50
Power [W]

Power [W]
30 30
40 40
20 20
Power [W]

Power [W]
30 30
10 10
20 20
0 0
100 5 10 15 20 25 100 5 10 15 20 25
Voltage [V] Voltage [V]
0 0
0 9: Characteristic
Figure 5 10 I-V 15
and P-V20curves 25 0 irradiation
for (a) varying 5 10 and 15
levels T = 25◦ C,
20 (b) varying
25
Voltage [V] 2 Voltage [V]
temperature levels and S = 1000W/m .

Derive expressions for enl and φ. Note, this model form is not arbitrary nor eye-balled from data.
(You are too skilled to accept such amateur approaches). It comes from an equivalent circuit
model detailed in Chapter 5 of Masters [6].

8 Offline Techniques

Offline model identification techniques consider scenarios where data is collected from the energy
system as a batch set. This data is then used for model identification, separate from the dynamic
evolution of the energy system itself. As a result, we do not recursively estimate the parameter
values as data arrives. Instead, we fit the model parameters in a single computational procedure.
In this section, we define this problem category and a general formulation. However, we do not
provide any specific algorithms and refer the reader to the optimization literature.
Consider the nonlinear-in-the-parameters model

y = f (u; θ) (78)

where y ∈ Rny is the model output, θ ∈ Rnθ is the true parameter vector, and u ∈ Rnu is a

Revised March 30, 2018 | NOT FOR DISTRIBUTION Page 22

CE 295 — Energy Systems and Control Professor Scott Moura — University of California, Berkeley

measured input vector that may or may not be present in the model. We can then formulate an
estimated output ŷ from estimated model parameters θ̂ as follows

ŷ = f (u; θ̂). (79)

With this notation, we can formulate the generic offline model identification problem as an opti-
mization task

minθ̂ ky − ŷk = ky − f (u; θ̂)k (80)

subject to θ̂ ∈ S (81)

where k · k in (80) is an appropriate norm (e.g. sum of squares, maximum magnitude, integrated
over time) and S in (81) is a set of feasible values for the unknown parameters θ̂.
Depending on the mathematical structure and size of the cost function and constraints (80)-
(81), various optimization techniques can be pursued. We provide two categories of techniques
and associated keywords that interested readers can search for further information

• [Gradient-based Methods] include the classical optimization techniques, such as gradient-

descent, quasi-Newton methods, and convex programming. These methods generally con-
tain attractive convergence and optimality properties. However, they require the cost/con-
straints to be differentiable in the variable θ̂.

• [Non-Gradient-based Methods] are pursued when it is intractable to compute (analytically

or numerically) the gradient. Several common methods include Genetic Algorithms, Particle
Swarm Optimization, and Stochastic Optimization. These algorithms typically require less
assumptions about the model structure, yet sacrifice computational efficiency and conver-
gence/optimality properties.

9 Sensitivity Analysis

Sensitivity analysis addresses the identifiability question. That is, this analysis is a structured way
to determine which parameters are identifiable. Moreover, it quantifies the level of identifiability,
allowing one to rank which parameters have the greatest influence on the output of interest. This
analysis is specifically oriented toward nonlinear-in-the-parameters models, where the output does
NOT have a linear dependence on the parameters. Despite our previous presentation of gradient,
least squares, and nonlinear least squares algorithms, sensitivity analysis is a necessary first step
in nonlinear parameterized models of energy systems.
Consider the nonlinear-in-the-parameters model

y = f (u; θ) (82)

Revised March 30, 2018 | NOT FOR DISTRIBUTION Page 23

CE 295 — Energy Systems and Control Professor Scott Moura — University of California, Berkeley

where y ∈ Rny is the model output, θ ∈ Rnθ is the parameter vector, and u ∈ Rnu is a measured
input vector that may or may not be present in the model. The sensitivity derivative of the model
outputs to the parameter vector is defined as

∂f
S = (u; θ) ∈ Rny ×nθ , (83)
∂θ  
∂f ∂f ∂f
S = [S1 , S2 , · · · , Snθ ] = (u; θ), (u; θ), · · · , (u; θ) (84)
∂θ1 ∂θ2 ∂θnθ

Each column of the sensitivity derivative, which we denote Si represents one parameter’s sen-
sitivity in all outputs, and is viewed as the sensitivity direction for the corresponding parameter.
Since the model is nonlinear in the parameters, global identifiability cannot be proven. That is, the
parameter’s sensitivity in an output may depend on the particular parameter values. Moreover, it
may also depend on the input vector u. In addition, only a subset of the parameters may be iden-
tifiable. This occurs when linear dependence exists between columns of the sensitivity derivative.
Intuitively, this means an output produces nearly identical reactions when two different parameters
are perturbed. As a result, it’s essentially impossible to differentiate between these parameters

y2 y2 y2

θ3 θ3 θ3
θ1 θ1 θ1

y1 y1 y1

θ2 θ2 θ2

Figure 10: Graphical description of Linear Dependence in Parameter Sensitivity

by observing the output. These complications motivate the need for an efficient and structured
method, which performs the following tasks:

1. Rank parameter sensitivity: This allows us to focus on identifying the most influential pa-
rameters.

2. Identify linear dependence between parameter sensitivities: This allows us to under-

stand which combinations of parameters are not identifiable, and focus only on the identifi-
able subset.

For the purposes of presenting this methodology, consider nθ = 4 parameters and an arbitrary
number of output data ny . Then the sensitivity derivative is composed of S = [S1 , S2 , S3 , S4 ] ∈

Revised March 30, 2018 | NOT FOR DISTRIBUTION Page 24

CE 295 — Energy Systems and Control Professor Scott Moura — University of California, Berkeley

Rny ×4 . A particular decomposition of matrix S T S reveals both sensitivity ranking and linear de-
pendence. Let S T S = DT CD where
   hS1 ,S2 i hS1 ,S3 i hS1 ,S4 i 
kS1 k 0 0 0 1 kS1 kkS2 k kS1 kkS3 k kS1 kkS4 k
   hS2 ,S1 i hS2 ,S3 i hS2 ,S4 i 
 0 kS2 k 0 0   kS2 kkS1 k 1 kS2 kkS3 k kS2 kkS4 k

D= , C=  , (85)
hS3 ,S1 i hS3 ,S2 i hS3 ,S4 i

 0 0 kS3 k 0 


 kS3 kkS1 k kS3 kkS2 k 1 kS3 kkS4 k


hS4 ,S1 i hS4 ,S2 i hS4 ,S3 i
0 0 0 kS4 k kS4 kkS1 k kS4 kkS2 k kS4 kkS3 k 1
√
where k · k denotes the Euclidian norm (kxk = xT x) and h·, ·i is the inner product (hx, yi =
xT y). This decomposition can be generated by performing Gram-Schmidt orthonormalization of
the sensitivity column vectors. This process coincides with the economical QR decomposition with
permutation (see qr command in Matlab). Now we discuss how to interpret this decomposition.

1. Rank parameter sensitivity: Diagonal matrix D provides a direct quantification of parame-

ter sensitivity. Sorting the elements kSi k in descending order provides the desired sensitivity
ranking. Note that one should normalize the parameters and outputs across a suitable range
to ensure a fair comparison.
hSi ,Sj i
2. Linear Dependence: By the Cauchy Schwarz inequality −1 ≤ kSi kkSj k ≤ 1. This has the
hSi ,Sj i
interpretation that values of kSi kkSj k near -1 or 1 imply strong linear dependence between
parameters θi and θj , whereas values near zero imply orthogonality. In other words, when
the off-diagonal element of matrix C is near -1 or 1, then the corresponding pair of param-
eters exhibit strong linear dependence, and are difficult to identify separately. When the
off-diagonal element is near zero, then the corresponding pair of parameters exhibit strong
linear independence, and can be easily identified separately. This concept is depicted graph-
ically in Fig. 10. The sensitivity vectors corresponding to θ1 and θ2 have an inner product
near -1, implying strong linear dependence. In contrast, the sensitivity vectors corresponding
to θ1 and θ3 have an inner product near 0, implying strong linear independence. Thus, θ1
and θ2 are difficult to identify separately, whereas θ1 and θ3 are easier to identify separately.

Example 9.1 (State-of-Health Parameter Estimation in Batteries). This example considers the
identification of uncertain parameters in a nonlinear battery model based upon electrochemistry
[7]. The nonlinear in the parameters model has the form:
   
−1 θ2 −1 θ3
V (t) = OCV + sinh √ I(t) − sinh √ I(t) + θ4 I(t) (86)
θ1 θ1

where V (t) is battery voltage, OCV is a known open circuit potential, I(t) is battery current, and
θ1 , θ2 , θ3 , θ4 are unknown electrochemical model parameters. Interestingly, identification of these
unknown parameters corresponds directly to assessing the battery’s age - specifically capacity
fade and power fade. For more details consult [7].

Revised March 30, 2018 | NOT FOR DISTRIBUTION Page 25

CE 295 — Energy Systems and Control Professor Scott Moura — University of California, Berkeley

An example of matrices D and C is provided in (87). This example analyzes parameter sensi-
tivity for a typical electric drive vehicle charge/discharge cycle.
   
2.827 0 0 0 1 −0.3000 0.2908 0.2956
   
 0 4.501 0 0   −0.3000 1 −0.9801 −0.9805 
D= , C=
 .

 0 0 3.958 0 
  0.2908 −0.9801 1 0.9322 

0 0 0 15.103 0.2956 −0.9805 0.9322 1
(87)
First, it is evident from D that the parameter sensitivity ranking is as follows: θ4 , θ2 , θ3 , θ1 . Accu-
rate identification of these parameters should be prioritized accordingly. In addition, strong linear
dependence exists between θ2 , θ3 , θ4 . This means perturbations in any of these three parameters
produces almost identical changes in our measured data. Consequently, attempting to individually
separate these variables is futile. An infinite number of linear combinations of these parameters
may satisfy our data. In this case, we should only identify the most sensitive parameter among
these linearly dependent subset. As a result, we conclude the model is partially identifiable and
focus on identifying the subset θ1 , θ4 only.
In this particular application, it turns out θ1 represents - exactly - charge capacity in the battery
and θ4 represents - exactly - internal resistance. Monitoring changes in these parameters yields
directly quantifications of charge and power fade.

Exercise 7 (Photovoltaic Array Learning - continued). blah blah...

10 Notes

Model identification is a rich and deep area, in both theory and applications. In this chapter, we
attempt to provide students with an introduction, several popular tools, and motivation for fur-
ther study. The quintessential reference for system identification is Lennart Ljung’s textbook [1].
Model identification is also commonly used in the context of control, in what is generally known as
“adaptive control.” Excellent textbooks in this area include Ioannou & Sun [2] and Krstic, Kanel-
lakopoulos, & Kokotovic [8]. Finally, Lund & Foss give a particularly excellent article on sensitivity
analysis in nonlinear parameterized models [9].

References
[1] L. Ljung, System identification. Wiley Online Library, 1999.
[2] P. Ioannou and J. Sun, Robust Adaptive Control. Prentice-Hall, 1996.
[3] J. Pukrushpan, H. Peng, and A. G. Stefanopoulou, “Control-oriented modeling and analysis for auto-
motive fuel cell systems,” Journal of dynamic systems, measurement, and control, vol. 126, pp. 14–25,
2004.

Revised March 30, 2018 | NOT FOR DISTRIBUTION Page 26

CE 295 — Energy Systems and Control Professor Scott Moura — University of California, Berkeley

[4] J. Pukrushpan, A. Stefanopoulou, and H. Peng, “Control of fuel cell breathing,” Control Systems, IEEE,
vol. 24, no. 2, pp. 30–46, Apr 2004.

[5] S. J. Moura and Y. A. Chang, “Lyapunov-based Switched Extremum Seeking for Photovoltaic Power
Maximization,” Control Engineering Practice, vol. 21, no. 7, pp. 971 – 980, 2013.

[6] G. M. Masters, Renewable and efficient electric power systems. John Wiley & Sons, 2013.

[7] S. J. Moura, N. Chaturvedi, and M. Krstic, “Adaptive PDE Observer for Battery SOC/SOH Estimation via
an Electrochemical Model,” ASME Journal of Dynamic Systems, Measurement, and Control, vol. 136,
no. 1, pp. 011 015–011 026, Oct 2014.

[8] M. Krstic, P. V. Kokotovic, and I. Kanellakopoulos, Nonlinear and adaptive control design. John Wiley
& Sons New York, 1995.

[9] B. F. Lund and B. A. Foss, “Parameter ranking by orthogonalization - Applied to nonlinear mechanistic
models,” Automatica, vol. 44, no. 1, pp. 278 – 281, 2008.

11 Appendix: Summary of Adaptive Laws

In this section, we summarize various model identification algorithms which are particularly useful in prac-
tice. These algorithms fall into two categories: Gradient Algorithms (Table 2) and Least Squares Algorithms
(Table 3).

Adaptive laws with Normalization: Throughout these tables, we consider adaptive laws with normaliza-
tion. These laws are identical to their non-normalized variants, except the estimation model is altered to
 = (z − ẑ)/m2 , where m2 is known as the normalization signal. This endows the algorithms with additional
boundedness properties, and generally provides better numerical conditioning.

Least-Squares with Forgetting Factor: Sometimes it is desirable to discount past data, relative to more
recent data. This is particularly useful for model parameters that might slowly change over time. The least
squares with forgetting factor is designed under this paradigm, where β is known as the forgetting factor.
When β = 0, then one recovers the pure least squares method. As β increases towards ∞, then one
approaches an instantaneous cost function.

Adaptive Laws with Projection: The parameter estimates θ̂ are generated by solving recursive opti-
mization problems, with appropriately formulated cost functions. In practice, one may wish to constrain
parameter estimates within a physically meaningful set. For example, heat transfer coefficients cannot
be negative, meaningful thermal resistances have upper and lower limits. Motivated by this situation, we
can add a convex constraint set to the optimization formulation. Mathematically, we enforce the constraint
θ̂ ∈ S = {θ̂ ∈ Rn | g(θ̂) ≤ 0}, where set S represents our convex constraint set that is defined by a function
g(·) that maps parameter estimate vectors to scalars. In this algorithm, we additionally require that g(·) is
a so-called “smooth” function. This means that g(·) and all its derivatives are continuous functions. This
is a mathematical technicality, but can be easily achieved in practice. A common selection for the convex
constraint set S is an ellipsoid, which is indeed parameterized by a smooth function g(·). In the algorithm
listed in Table 3, int(S) signifies the interior of set S, and δ(S) signifies the boundary of set S.

Revised March 30, 2018 | NOT FOR DISTRIBUTION Page 27

CE 295 — Energy Systems and Control Professor Scott Moura — University of California, Berkeley

Table 2: Gradient Algorithms

Parametric Model z = θT φ

Estimation model ẑ = θ̂T φ

z−ẑ
Normalized Estimation model = m2

Based on Instantaneous Cost

Cost J(θ̂) = 12 2
˙
Adaptive Law θ̂ = Γφ

Design variables m2 = 1 + αφT φ, α > 0, Γ = ΓT  0

Based on Integral Cost

Rt
Cost J(θ̂) = 1
2 0 e−β(t−τ ) 2 (t, τ )m2 (τ )dτ
˙
θ̂ = −Γ(Rθ̂ + Q)

Adaptive Law φφT

Ṙ = −βR + m2
, R(0) = 0

zφ
Q̇ = −βQ + m2
, Q(0) = 0

Design variables m2 = 1 + αφT φ, α > 0, β > 0, Γ = ΓT  0

Gradient Law with Projection

Cost J(θ̂) = 12 2

subject to θ̂ ∈ S = {θ̂ ∈ Rn | g(θ̂) ≤ 0, g : Rn → R is smooth fcn}






 Γφ if θ̂ ∈ int(S), OR

˙

Adaptive Law θ̂ = Proj(Γφ) = if θ̂ ∈ δ(S) and (Γφ)T ∇g ≤ 0



 T
 Γφ − Γ ∇g∇g

Γφ otherwise

∇g T Γ∇g

Design variables m2 = 1 + αφT φ, α > 0, Γ = ΓT  0, g(·) smooth

Revised March 30, 2018 | NOT FOR DISTRIBUTION Page 28

CE 295 — Energy Systems and Control Professor Scott Moura — University of California, Berkeley

Table 3: Least Squares Algorithms

Parametric Model z = θT φ

Estimation model ẑ = θ̂T φ

z−ẑ
Normalized Estimation model = m2

Pure Least Squares with Normalization

Rt T (t)φ(τ )]2
Cost J(θ̂) = 12 0 [z(τ )−mθ̂ 2 (τ )
dτ
˙
θ̂ = P φ
Adaptive Law
T
Ṗ = −P φφ
m2
P, P (0) = P0

Design variables m2 = 1 + αφT φ, α > 0, P0 = P0T  0

Least Squares with Forgetting Factor

Rt T (t)φ(τ )]2
h iT h i
Cost J(θ̂) = 12 0 e−β(t−τ ) [z(τ )−mθ̂ 2 (τ )
dτ + 1 −βt
2 e θ̂(t) − θ̂ 0 Q0 θ̂(t) − θ̂ 0

˙
θ̂ = P φ
Adaptive Law
T
Ṗ = βP − P φφ
m2
P, P (0) = P0 = Q−1
0

Design variables m2 = 1 + αφT φ, α > 0, β > 0, P0 = P0T  0

Least Squares with Forgetting Factor and Projection

Rt T (t)φ(τ )]2
h iT h i
Cost J(θ̂) = 12 0 e−β(t−τ ) [z(τ )−mθ̂ 2 (τ )
dτ + 1 −βt
2 e θ̂(t) − θ̂ 0 Q0 θ̂(t) − θ̂ 0

subject to θ̂ ∈ S = {θ̂ ∈ Rn | g(θ̂) ≤ 0, g : Rn → R is smooth fcn}


P φ if θ̂ ∈ int(S), OR





˙

Adaptive Law θ̂ = Proj(P φ) = if θ̂ ∈ δ(S) and (P φ)T ∇g ≤ 0


 T
 P φ − P ∇g∇g

T P ∇g P φ otherwise

 ∇g
T
βP − P φφ 2 P if θ̂ ∈ int(S), OR


m




Ṗ = if θ̂ ∈ δ(S) and (P φ)T ∇g ≤ 0




 0 otherwise


Design variables m2 = 1 + αφT φ, α > 0, β > 0, P0 = P0T  0, g(·) smooth

Revised March 30, 2018 | NOT FOR DISTRIBUTION Page 29