Note: Descriptions are shown in the official language in which they were submitted.
CA 02620743 2008-01-10
-1-
METHODS FOR OPTIMIZING SYSTEM MODELS AND THE LIKE USING
LEAST ABSOLUTE VALUE ESTIMATIONS
Technical Field
The invention disclosed herein relates to methods for fitting curves to sets
of data points
(e.g. experimentally determined data points). Particular embodiments of the
invention
provide methods for optimizing solutions to least absolute value (LAV) curve
fitting
problems.
Background
The scientific method depends on conducting experiments and analyzing the data
obtained from such experiments. Often, the data analysis involves statistical
methods.
When analyzing experimentally obtained data, it is often helpful to fit the
experimentally
obtained data to a curve. This process may be referred to as "curve fitting".
Curve fitting
typically involves defining a "curve" (i.e. a function sometimes referred to
as a "model")
that relates one data value to one or more other data values and selecting
parameters of
the curve such that the curve estimates the relationship between the data
values. By way
of non-limiting example, selection of the parameters of the curve may involve
selection
of coefficients of a polynomial function.
The curve established in the curve fitting process is typically an estimate,
and not an
exact representation, of the relationship between the known data values. The
curve fitting
process may involve the use of a cost function which provides a metric as to
how closely
the curve determined in the curve fitting process estimates the relationship
between the
data values. Curve fitting processes often seek to determine a curve which
minimizes (or
maximizes) the cost function. Minimization (or maximization) of the cost
function may
be referred to as "optimizing" the fitted curve.
One well known example of a cost function is the "least squares" cost
function. The least
squares cost function represents a sum of the squares of the residuals - i.e.
a sum of the
squares of the differences between values generated by the fitted curve and
correspondirig
CA 02620743 2008-01-10
-2-
values from the data. Optimizing a fitted curve using a least squares cost
function is a
well known problem for which there are a number of known solutions. Examples
of such
solutions include the so called "Newton's method" and "gradient descent
method".
There are many drawbacks with the least squares cost function and the
associated
optimization methods. For example, iterative optimization is generally
required to
minimize the cost functions associated with non-linear models. Iterative
optimization is
typically performed using suitably programmed computer(s) or similar data
processor(s).
Using a least square cost function and the gradient descent method requires a
relatively
1 o large number of iterations for the optimization to converge. Optimization
techniques
involving a large number of iterations are computationally expensive (i.e.
requiring
relatively large amounts of processing time and/or computational storage).
Another
drawback associated with least squares optimization processes is that they are
relatively
non-robust to outlying data values. Accordingly, if there is an obvious error
with one or
more of the data values (e.g. because a number was mistyped into a keyboard or
a
measurement apparatus was uncalibrated at the time of measurement), the
resultant
model obtained using least squares optimization could be skewed.
There is a desire to provide optimization methods for curve fitting which are
relatively
2o easy to perform and which ameliorate some of the drawbacks associated with
prior art
techniques.
Brief Description of the Drawings
In drawings which schematically illustrate non-limiting embodiments of the
invention:
Figure I is a depiction of an exemplary environment suitable for practicing
various embodiments of the invention;
Figure 2 shows a method for determining a model of the Figure 1 system under
consideration according to a particular embodiment of the invention;
Figure 3 shows a method for optimizing a curve/model to a set of data values
according to a particular embodiment of the invention where the curve/model is
unconstrained;
CA 02620743 2008-01-10
-3-
Figure 4 shows a particular method for obtaining a values suitable for use
with
the optimization methods described herein;
Figure 5 shows a method for optimizing a curve/model to a set of data values
according to a particular embodiment of the invention where the curve/model is
constrained by particular constraint criteria;
Figure 6 shows a method for optimizing a non-linear curve/model to a set of
data
values according to a particular embodiment of the invention; and
Figure 7 shows a method for optimizing a non-linear curve/model to a set of
data
values according to another embodiment of the invention.
Description
One aspect of the invention provides a method for fitting a curve to a set of
data values
using a least absolute value (LAV) cost function. The set of data values may
comprise m
sets of data values. The method takes advantage of contraction mapping to
determine a
number n<m of individual equations which are interpolated by the curve to be
fitted,
where n corresponds to the number of parameters x, , x2,..., xõ to be
ascertained for the
curve to be fitted. The method then involves solving the n selected equations
to determine
the n parameters xõ xz,..., xõ of the curve to be fitted. Selection of these
parameters
ensures that the curve to be fitted minimizes the LAV cost function.
The methods may be used for linear and non-linear models with or without
constraints on
the model as described in more particular detail below.
Figure 1 schematically depicts an exemplary environment 10 in which the
particular
embodiments of the invention may be employed. Environment 10 comprises a
system of
interest 12. System 12 is provided with one or more inputs 14 and an output
16. In the
illustrated embodiment, system 10 comprises one output 16 and as many as n
inputs 14.
Output 16 may be measured by measurement apparatus 18, which, by way of non-
limiting example, may comprise one or more sensors, transducers or the like.
In other
3o embodiments, system 12 may comprise a plurality of outputs 16. In some
embodiments,
output(s) 16 are provided by system 12 in a form which can be determined
without
CA 02620743 2008-01-10
-4-
explicit measurement, in which case measurement apparatus 18 is not required.
In the
particular exemplary environment 10, it is desired to obtain a model 26 which
approximates the response of system 12 to inputs 14 - i.e. model predicts
output 16 for a
given set of inputs 14. Mode126 may also be referred to herein as a"curve" or
the like.
In the illustrated embodiment, model 26 is determined by processor(s) 20 which
perform(s) a methods according to particular aspects of the invention as
described in
more particular detail below. Processor(s) 20 may generally include any data
processing
device(s) capable of performing the functions described herein. By way of
example,
processor(s) 20 may comprise one or more suitably programmed computer(s), one
or
more suitably programmed microprocessors, one or more suitably programmed PLC
devices or the like and may include a combination of such data processing
devices.
In the illustrated embodiment, processor(s) 20 are configured to optionally
control inputs
14 via signal(s) 22 and processor(s) 20 are configured to receive the measured
value of
output 16 from measurement apparatus 18 via signal(s) 24. Suitable hardware
and/or
software (not explicitly shown) may be used to effect the communication of
signals 22,
24. In some embodiments, processor(s) 20 may be provided with information
relating to
inputs 14 and output 16 using other techniques. For example, inputs 14 may be
manually
controlled or controlled by a separate input controller and output 16 and/or
measurement
apparatus 18 may be manually read. The values of inputs 14 and output 18 may
be
provided to processor(s) 20 by some other suitable input means. By way of non-
limiting
example, the values of inputs 14 and output 16 may be keyboarded into
processor(s) 20
or read into processor(s) 20 from a disc drive.
Figure 2 schematically depicts a method 50 according to a particular
embodiment of the
invention which may be used in environment 10 to determine a system model 26.
Method
50 may be performed, at least in part, by processor(s) 20. Method 50 starts in
block 52
which involves obtaining the values of inputs 14. Inputs 14 may be provided to
processor(s) 20 using any suitable technique including, without limitation,
those
discussed above. Method 50 then proceeds to block 54 which involves obtaining
the
CA 02620743 2008-01-10
-5-
values of output 16 corresponding to each set of n inputs. Outputs 16 may be
provided to
processor(s) 20 using any suitable technique including, without limitation,
those
discussed above. In total, blocks 52 and 54 may involve obtaining m sets of
inputs 14 and
m corresponding outputs 16.
Method 50 then proceeds to block 56 which involves obtaining a form of model
to be
considered. The form of model may be provided to processor(s) 20 in any
suitable
manner. By way of non-limiting example, the form of model may be provided to
processor(s) 20 via software instructions. Such software instructions may be
provided to
1o processor(s) 20 by any suitable input means and may provide processor(s) 20
with the
form of model when the software instructions are executed by processor(s) 20.
In another
non-limiting example, processor(s) 20 may be provided with a user interface
that queries
a user to input the form of model in block 56. In the particular embodiment of
method 50
illustrated in Figure 2, system model 26 is determined by a number n of model
parameters.
Method 50 then proceeds to optional block 58 which involves obtaining any
constraint(s)
which may need to be considered. Optional constraints may be provided to
processor(s)
in any suitable manner. By way of non-limiting example, the optional
constraints may
2o be provided to processor(s) 20 via software instructions. Such software
instructions may
be provided to processor(s) 20 by any suitable input means and may provide
processor(s)
20 with the constraints and/or the form of constraint when the software
instructions are
executed by processor(s) 20. In another non-limiting example, processor(s) 20
may be
provided with a user interface that queries a user to input any optional
constraints in
block 58.
Method 50 then proceeds to block 60 which involves determining the model
parameters
using LAV techniques. It will be appreciated that in the illustrated
embodiment of
method 50, detenmining the model parameters is sufficient to specify model 26
of system
12. Particular embodiments illustrating a number of techniques suitable of use
in block 60
CA 02620743 2008-01-10
-6-
for detenmining model parameters using LAV techniques are described in more
detail
below.
Method 50 then proceeds to optional block 62 which involves outputting the
model
parameters determined in block 60. It will be appreciated that in the
illustrated
embodiment of method 50, outputting the model parameters is sufficient to
specify model
26 of system 12.
The Linear Measurement Equation
We define a set of data values which include:
(i) a data matrix H (where H E R 's" is of maximal rank); and
(ii) a mxl data vector z,
where m>n. The elements of data matrix H and the elements of the data vector z
may
represent experimentally obtained data, for example. In this description, we
refer to the
individual elements of the data matrix H as ay, where i represents the row
andj represents
the column of the element au. In some embodiments, the elements of each row of
data
matrix H (e.g. a a,z, ... aiõ) may represent "input" or "controllable" or
"independent" or
,~,
"explanatory" variables (e.g. inputs 14 of environment 10 (Figure 1)) and the
element in
each row of data vector z (e.g. z,) may represent a corresponding "output" or
"measured"
or "dependent" or "response" variable (e.g. output 16 of environment 10
(Figure 1)).
Despite the use of this language to describe the data variables H and z, in
general,
causality between data matrix H and data vector z is not a required feature of
the
invention disclosed herein.
In accordance with particular embodiments of the invention, we assume that
each
element z; of output data vector z is a linear combination of n state
variables xi, x2, ... xõ
to be estimated corrupted by a corresponding disturbance v;. The n state
variables xl,
xZ, ...xõ may be represented by an nx I state variable vector x and the m
disturbances v;
may be represented by an uncertain mx 1 disturbance vector v. The state
variable vector x
CA 02620743 2008-01-10
-7-
may represent mode126 (Figure 1) of system 12 and the n state variables xl,
x1, ...xõ may
represent the parameters of mode126.
We may then express the relationship between z. H, x, and v in the following
form
z=Hx+v (1)
For the sake of convenience, equation (1) may be referred to as a "measurement
equation" and each of the individual elements z; of the output data vector z
may be
referred to as "measurements".
The i`h measurement z; in measurement equation (1) can be written in component
form as:
z, = H, x+ v, i=1,2,...,m (2)
where H; E R'x" and x is an n-dimensional state vector made up of parameters
xf, xz, ... x,,. In accordance with particular embodiments of the invention,
it is desired to
estimate the parameters xi, x2, ...xõ which determine the n-dimensional state
vector x(e.g.
model 26) using output data values z, (e.g. from outputl6) and input data
values H,
(e.g. from inputs 14). Method 50 (Figure 2) represents a particular embodiment
of a
method for obtaining such an estimate.
Least Absolute Value Estimator
One aspect of the invention provides new methods which make use of least
absolute
value (LAV) residuals to estimate the n-dimensional state vector x (i.e. the
parameters
xl, x2, ... xõ) in the measurement equation (1). By way of non-limiting
example, such
methods may be used to implement block 60 of method 50 (Figure 2). As used
herein, an
individual LAV residual represents the absolute value of the difference
between the value
predicted by the proposed (or estimated) function (e.g. mode126) and the
corresponding
measurement z;. In particular embodiments of the invention, methods are used
to
determine the value of x (i.e. the parameters xi, X2, ...xõ) which minimizes a
cost function
representing the sum of the LAV residuals over the m individual measurements
(zi, Z2,
...zõ)
CA 02620743 2008-01-10
-8-
It may be shown that:
Theorem if the column rank of the matrix H E R mrn in (1) is k<_ n < m, then
there exists
a best LAV approximation (referred to herein as HxL,,,, ) which interpolates
at least k
points of the set z=[z1, zZ , . . . , Z. I T in (1).
In other words, if we have m measurements zl, zz, ... zm (e.g. m outputs 16)
and input data
values (e.g. from inputs 14) corresponding to each of the m measurements to
form
H E Rmzn and if the rank of H is k< n < m, then the best LAV approximation Hx,-
,,,.
interpolates at least k of the m measurements zl, zz, ... zm. A particular
measurement z; is
l0 interpolated when its corresponding disturbance v; is zero in measurement
equation (1).
Accordingly, if the conditions of the above theorem are met, then the
disturbance vector v
of measurement equation (1) includes k elements for which v;=0. This
interpolation of the
best LAV approximation Hx,,,. contrasts with the least squares approximation
which
does not necessarily pass through any of the individual measurements zi, Z2,
... z,õ.
For the measurement equation (1), the LAV cost function to be minimized may be
written as
m m
J(x)=jIv,l =Erow,lz-Hxl (3)
;=1 ,=1
In an n-dimensional linear space, the so called p - norm of a vector v [
= v, , vZ, ... , v]~_ õ
can be defined by means of the formula
n n
11vII p - ( I lvk I P < 00 p =1,2,... (4)
k=1
For p=1 and k=1,..., n equation (4) can be written as
11vII =Iv,I +IvZ1 +...+Ivn 1 (5)
where II=II represents the I-norm in this description.
CA 02620743 2008-01-10
-9-
On the basis of the norm axioms, we know Ilv_II = 0 if and only if v = 0.
Consequently, the
cost function in equation (3) is a minimum when
v, =(z, - H;x) = 0 for i= 1,2,...,m (6)
Equation (6) characterizes an over-determined system of equations, since m>n.
Optimization using LAV estimation is advantageous over prior art techniques
because, as
pointed out by the above-presented Theorem, a number k of the individual
measurements
are interpolated by the approximation Hx,,,,. As described in more detail
below, this
feature results in outlying measurements being rejected by the optimization
methods of
1o the invention. The rejection of outlying measurements can be seen by
realizing that the
residual error due to an outlying measurement will be significant.
Consequently, the
outlying measurement will not be one of the k measurements interpolated by the
approximation Hx,_4y. .
Contraction Mapping
The usual representation for any optimization problem is a correspondence
between the
measurement vector z and the estimated state vector x. The measurement vector
z and
the estimated state vector x are identified as belonging to a metric space,
endowed with a
mathematical structure suited to the type of analysis desired. A metric space
X which has
the property that every Cauchy sequence in X has a limit in X is said to be
"complete".
Furthermore, if a normed linear space is complete, it is called a "Banach"
space.
In particular embodiments of the invention, the space best suited for our
analysis is a
Banach space. As described in more detail below, a Banach space is convenient,
because
it facilitates the application of contraction mapping when we desire to
determine the
optimal value of x in (6).
Consider the set of ordered n-tuples of real numbers with the distance
function
n
d(x, y) = E lxk - yk 1 (7)
k=1
CA 02620743 2008-01-10
-10-
Ix-YII
which forms the metric space 1, and corresponds to the norm given by equation
(4) with
p=1. It is readily seen that the function Ilxll, defined on R" by (7) forms a
norm on the
Banach space R".
Consider the mapping f: R" -> R", x H Ax given by the system of linear
equations
n
x; = Ea,,xj +b; i=1,2,...,n (8)
J=I
where a, are the elements of the matrix A and compare this system of equations
to the
measurement equation (1). A mapping A, of R", into itself is referred to as a
"contraction mapping" if there exists a number a<1 such that
d(Ax, Ay) S ad (x, y) t>x, y E R" (9)
If a matrix A is a contraction mapping, then we can apply a method of
successive
approximations to determine the solution of the equation
x=Ax (10)
It can be shown that every contraction mapping defined in the Banach space
(R",II=II) has
one and only one fixed point, i.e., the equation Ax; = x;+, has one and only
one solution
x =1imn~m xn =0 (11)
Furthermore, it can be shown that the distance between two values of x, is
given by
d(xn~x,n) :~ a"d(XO,x, )l(1-a) (12)
where a<l. It may be concluded, therefore, that the sequence forms a
Cauchy sequence since the right hand side of (12) tends to zero as n tends to
infinity.
When equation (12) is true, then we can always approximate the solution of
equation (11)
as indicated in (13). We can limit the solution of (13) by excluding only the
last few
CA 02620743 2008-01-10
-11-
consecutive terms in (13). However, the larger the number of terms we use in
(13) the
more accurate is the solution given by (13).
The solution of (10), obtained by successive approximations, is given by
x = xa +(x, -xo)+...+(xõ -xn_,) (13)
and inspection of (12) clearly indicates that (13) converges rapidly if a 1.
Having
a 1 is desirable, since we want to perfonn the smallest number of iterations
possible,
in order to obtain a good estimate of the solution x.
We now apply these contraction mapping results to the system of equations
given by (8)
using the distance function (7). It can then be shown (see Appendix A) that
(8) is a
contraction mapping if
1 l a;~ l<_ a< 1 `dj E{1,..., n} (14)
that is, only if (8) is seen as a mapping of the form
f:(R ,11'1)`+ (R",Il'll),xHAx+b
Particular Embodiment for Unconstrained Optimization
Figure 3 illustrates a particular embodiment of a method 100 for optimizing a
curve/model x defined by n parameters x,, x2,..., x,, for a set of data values
H, z according
to a particular embodiment of the invention. For example, method 100 may
involve
estimating n parameters x, , xz,..., xõ which define a mode126 of system 12
(Figure 1).
Method 100, which comprises determining the n model parameters x, , xz ,...,
xõ using
LAV techniques, may be performed in block 60 of method 50 (Figure 2).
It will be recalled from the above-discussed Theorem that if the rank of the
matrix H in
measurement equation (1) is k, the best LAV approximation Hx,,,,, interpolates
k of the
measurements zi, z2, ... zm in (1). For these k interpolated measurements z;,
the
corresponding v, = 0 will be zero in equation (3). As described in more detail
below,
method 100 makes use of the m measurements zi, z2, ... z,õ in (1), but
involves selection
CA 02620743 2008-01-10
-12-
of a desired subset of n measurements to be interpolated using relation (14)
where we
assume that k=n in the Theorem stated above.
Referring to Figure 3, block 102 involves obtaining data values corresponding
to data
matrix H and measurements z. In particular embodiments, the data values H, z
may
comprise experimentally determined data obtained from system of interest 12
(e.g. a mxn
input matrix H of inputs 14 and a mxl output vector of measured outputs 16).
The data
obtained in block 102 of Figure 3 may correspond to inputs 14 of system 12 and
corresponding outputs 16 of system 12 obtained in blocks 52, 54 of method 50
(Figure2).
System 12 may generally include any system from which data may be measured or
otherwise ascertained. By way of non-limiting example, system of interest 12
may be an
electrical circuit, a wireless communication system, a chemical system, a
mechanical
system, a financial system, a theoretical system or any other suitable system.
In some
embodiments, the data values H, z may comprise m data sets, where each data
set
comprises one or more independent (explanatory) variables (e.g. inputs 14) and
a
corresponding measured dependent (response) variable (e.g. outputs 16). By way
of non-
limiting example, in an electrical circuit system, a set of data values may
comprise
independent (explanatory) variables which include the amplitude of voltage
and/or
current inputs to the system and a measured dependent (response) variable that
is the
measured voltage output of the system. As another non-limiting example, in the
case of a
chemical system, a set of data values may comprise independent (explanatory)
variables
which include the temperature of the reaction chamber and/or the concentration
of a
reactant X and a measured dependent (response) variable which may include the
concentration of the reaction product Y.
While the block 102 data values in particular embodiments of the invention
include
experimentally determined data, this is not necessary. In general, the block
102 data
values may comprise m sets of data values and each set of data values may
comprise n+1
individual data elements. These data values may be obtained in any suitable
manner.
CA 02620743 2008-01-10
-13-
Method 100 then proceeds to block 104, which involves generating a measurement
equation of the form described above in (1). The data matrix H and the data
vector z are
constructed from the block 102 data values. The block 102 data values may
comprise m
sets data values and each such set may comprise n+1 individual data elements.
In block
104, each set of n+1 individual data elements from block 102 is used to form
one row H;
(consisting of n individual data elements (a; f, arz, ... a;")) of the data
matrix H and one
row z; (consisting of a single data element) of the data vector z. Since the
block 102 data
values comprise m sets of data, the data matrix Hof the block 104 measurement
equation
is an mxn dimensional matrix and the data vector z of the block 104
measurement
equation is a mxl dimensional vector.
The state variable vector x in the measurement equation (1) is a nxl vector
and the
individual elements xi, X2, ... x, of the state variable matrix x represent
the variable
parameters (as yet unknown) of the curve/model to be fitted (e.g. the
parameters (as yet
unknown) of model 26). In general the parameters xl, x2, ... x" may correspond
to any
suitable parameters of any suitable curve. By way of non-limiting example, the
curve to
be fitted may be a polynomial curve of the form f(b) = x, + x2b + x3b2 +...
x"b"-' such
that the individual elements xl, x1, ... x" of the state variable matrix x
represent the
coefficients of the polynomial. As another non-limiting example, the curve to
be fitted
may be an exponential curve of the form f(b) = x,e'=h . The disturbance vector
v is a mxl
vector whose individual elements vi, vz, ... v,õ represent the (as yet
unknown) error
between the fitted curve and the individual measurements zl, zz, ... z,õ.
It will be appreciated that the block 104 measurement equation (as represented
by (1))
actually comprises m individual equations having the form of equation (2).
Block 106 involves an optional inquiry into whether any of the m individual
equations in
the block 104 measurement equation (1) are identical to one another. If there
are no
identical equations (block 106 NO output), then method 100 proceeds to block
109 with
the full set of m individual equations. On the other hand, if there are one or
more identical
equations (block 120 YES output), then method 100 proceeds to block 108 where
such
CA 02620743 2008-01-10
-14-
identical equations are removed from the measurement equation to yield a
measurement
equation with a reduced set of mi individual equations (ml>n). Unless
explicitly
mentioned, it is assumed for the remainder of this description that the block
106 inquiry
is negative (i.e. that method 100 proceeds with a full set of m individual
equations). It
will be appreciated that this assumption yields no loss of generality,
provided that ml>n.
That is, the results obtained using a reduced set mj>n of equations are
obtained in a
similar manner to those obtained using a full set m of equations by replacing
m with mi.
Block 109 involves normalizing the block 104 measurement equation (1). The
block 109
normalization process involves dividing each of the m individual equations in
the block
104 measurement equation (1) by a suitable corresponding normalizing constant
such that
the data vector z becomes z=[z,, zZ,..., zm ]" =[0.1,0.1,...,0.1]" . To
achieve a data vector
z=[z,,zz,...,zm]' _[0.1,0.1,...,0.1]' the block 109 normalizing constants may
be
different as between the m individual equations in the block 104 measurement
equation
(1). It is a well known mathematical principal that dividing both sides of an
equation by a
constant does not alter the validity of the equation. The remainder of the
description of
method 100 assumes that we are referring to the normalized measurement
equation
obtained in block 109.
In the illustrated embodiment, block 110 comprises determining a set of
parameters
referred to as "a values" for the normalized H matrix. In general, for a
normalized
measurement equation having a mxn input matrix H. there are n a values for
each
possible combination of n rows of the normalized H matrix. For the purposes of
this
description, we adopt the following notation
aj.kd.... =laj,l+la, I +lar,I +... (14')
where:
= i corresponds to an a value calculated from the elements of the i'" column
of the
normalized H matrix;
= j, k, 1, ... corresponds to a row combination of the jr", kt",1`" . rows,
where the
number of rows in any particular row combination is n; and
CA 02620743 2008-01-10
-15-
ay corresponds to the element of the normalized H matrix at row x, column y.
In accordance with equation (14'), the a value a;ll..k Jis calculated by
summing the
...
absolute values of elements in the i`h column from thej'h, k`h, l'h ... rows
of the
normalized H matrix - i.e. (ap I+ lak, I+ la,; I+===. The total number of
elements a.,y of the
normalized H matrix used to calculate a particular a value (corresponding to
the number
of rows in any given row combination) is n and these n elements are always
selected from
the same column of the normalized H matrix. In some instances, an a value a
may
Y
be referred to in an abbreviated form a;, where the i refers to an a value
calculated from
the elements of the i'h column of the normalized H matrix and the particular
row
combination j, k, 1, ... used to calculate the a value is omitted from the
notation.
Figure 4 shows a method 110 for determining the particular a values of a given
normalized H matrix according to a particular embodiment of the invention. It
will be
appreciated that method 110 illustrated in Figure 4 represents one among many
suitable
methods which may be used to implement block 110 of method 100 (Figure 3) -
i.e. to
determine the a values of the normalized H matrix. Generally speaking, block
110 of
method 100 (Figure 3) may be implemented using any suitable method for
calculating
these a values. Method 110 starts in block 202, where the a, values are
calculated for all
of the possible non-ordered row combinations. of n rows which include row 1 of
the
normalized H matrix. That is, all combinations of a, I j=Oõ = a,jI+ lak~ I+
laõ I+- for
which j=1 and k,I,... E{2,...,mlk # 1# ...}.
Method 110 then proceeds to block 204 which involves calculating the az values
for all
of the possible non-ordered row combinations of n rows which include row I of
the
normalized H matrix. That is, all combinations of a Z I j kl = a,2+(ak2I+ Ia,Z
I+=== for
which j=1 and k,1,... E{2,...,mlk # 1:A ...}. After block 204, method 110
continues to
calculate the a values for the row I combinations in this manner, until it
calculates the aõ
values for the row I combinations in block 206. Block 206 involves calculating
CA 02620743 2008-01-10
-16-
aõ=la,õI+la,p,I+la,,,I+=== forwhichj=l and k,l,...e{2,...,mlk* 1~...}.Atthe
~. . .
conclusion of block 206, method 110 has calculated all of the a,, a2, ... aõ
(corresponding to columns 1, 2, ... n of the normalized H matrix) for all of
the possible
non-ordered combinations of n rows which include row 1(i.e. the row 1
combinations).
Method 110 then proceeds to block 208 which involves calculating the af values
for all
of the remaining possible non-ordered row combinations of n rows which include
row 2
of the normalized H matrix. The block 208 row combinations include row 2 but
do not
include row 1, as the a, values for all possible row combinations of n rows
that include
row I were calculated in block 202. The a, values calculated in block 208
include all
combinations of a, I j.k.1, = laj, I+ ak, I+ laõ I+=== for which j=2 and
k,1,... E~3,..., mlk :t- 1* ... }. Similarly, block 210 involves calculating
the a2 values for the
remaining possible row 2 combinations. The a2 values calculated in block 208
include all
combinations of a2 I~ k~ = laiZ I+ lak2 j + la,Z I+-=- for which j=2 and
k,1,... E13,..., mlk # 1~... }. Method 110 continues to calculate the a values
for the
remaining row 2 combinations until block 212, which involves calculating the
aõ values
for all combinations of aõI j,k.1, ., = l ajõ I+ l akõ I+ la,. I+ -- = for
which j=2 and
k,1,... E{3,..., ml k# 1# ... }. At the conclusion of block 212, method 110
has calculated (in
blocks 208, 210 ... 212) all of the a,, a2i ... aõ (corresponding to columns
1, 2, ... n of
the normalized H matrix)for all of the non-ordered combinations of n rows
which include
row 2 but do not include row 1.
It will be appreciated that the number of remaining row 2 combinations for
blocks 208,
210, 212 is fewer that the number of row I combinations in blocks 202, 204,
206. Method
110 continues in this manner until it reaches row m-n+l in block 214. When
method l 10
reaches row m-n+1, there is only one remaining combination of n rows for which
the a
values have not been calculated - i.e. the row combination m-n+1, m-n+2, ...
m.
Accordingly, block 214 involves the calculation of the a, value given by
CA 02620743 2008-01-10
- ] 7-
a,lj.k.1 =Ia,,I+Iak,I+laõI+--- forwhichj=m-n+l,k j+l, 1=k+1 ...,block216
involves the calculation of a2l j,k,1,.. = Ia,ZI +lak2l +lafzl +' =- for
whichj= m-n+l, k j+l,
l=k+ 1... and block 218 involves the calculation of a. l j.k,, = lajõ I+ lakõ
I+ la,j+=-- for
whichj= m-n+l, k j+l, 1=k+1 ....
At the conclusion of block 218, method 110 has calculated all of the a,, a2,
... an for all
of the possible non-ordered combinations of n rows. Method 110 is therefore
complete.
It will be appreciated that the number of a values calculated in block 110
(e.g. by method
110) is related to the values of m and n. It is well known in statistics that
the numbered to
non-ordered combinations of n elements (e.g. n rows) which may be selected
from a
group of m elements is given by m= m! . Therefore, in the illustrated
n n!(m-n)!
embodiment, block 110 comprises determining m! a, values, m! a2
n!(m-n)! n!(m-n)!
values ... 'n! aõ values for a total of nm! = m! a values.
n!(m-n)! n!(m-n)! (n-1)!(m-n)!
Returning to method 100 of Figure 3, after determining the a values in block
110, method
100 proceeds to block 115. Block 115 involves an inquiry into whether all of
the a values
determined in block 110 are less than unity. If any of the a values determined
in block
110 are greater than or equal to unity (block 115 NO output), then method 100
proceeds
to block 117, where the entire measurement equation is divided by a suitable
factor of 10
before returning to block l 10 to determine a new set of a values. In the
illustrated
embodiment, the loop of blocks 110, 115 and 117 continues until all of the m
values of a
are strictly less than unity. In this circumstance (block 115 YES output),
method 100
proceeds to block 130.
CA 02620743 2008-01-10
-18-
In some embodiments, the block 115 inquiry may be performed after the
calculation of
each a value (or at some other suitable interval prior to completing the
calculation of all
of the a values). If any one a value is greater than or equal to unity (block
115 NO
output), then the entire measurement equation is divided by a suitable factor
of 10 (block
117). Performing the block 115 inquiry after the calculation of each a value
(or at some
other suitable interval prior to completion of calculation of all of the a
values) may save
computational resources in cases where an a value greater than or equal to
unity is
determined early, since method 100 will not proceed to block 130 unless all of
the a
values are strictly less than unity and therefore there is no point in
calculating any further
lo a values using the current measurement equation if any one of the a values
is determined
to be greater than or equal to unity.
Block 130 involves selecting n individual equations from among the set of m
(or mi)
individual equations. The n individual equations are selected on the basis of
the row
n
combination for which a, = a, + a2 +...an is the lowest. The quantity
n
a, = a, + a Z+...aõ for a particular row combination may be referred as the a
sum for
;-~
that row combination. In the illustrated embodiment, block 130 involves
selecting the one
row combination with the lowest a sum and then selecting the n individual
equations
corresponding to the rows of the row combination with the lowest a sum. By way
of non-
limiting example, the block 130 selection of n individual equations may
involve summing
the a values obtained in block 110 (i.e. a, = a, + a 2+...a n) for each row
combination and then selecting the n individual equations corresponding to the
row
combination with the lowest a sum.
The result of the block 130 selection is a set of n individual equations. This
set of n
individual equations has 2n unknowns which include the corresponding
parameters
CA 02620743 2008-01-10
-19-
x,, xZ,..., x,, of the state variable vector x and the corresponding values
v,, v2,..., võ of the
disturbance vector v. However, on the basis of the Theorem described above, we
assume
that k=n and that the best LAV approximation Hx,,,,. interpolates the n
measurements
selected in block 130. Accordingly, for the n equations selected in block 130,
the
corresponding values of the disturbance vector v are all zero. We therefore
have reduced
the set of equations selected in block 130 to a set of n individual equations
with n
unknowns (i.e. the corresponding parameters x,, xZ,..., x,, of the state
variable matrix x).
The n equations selected in block 130 represent the n individual equations
with the
smallest contraction mapping constants (a values) in (9) and (14). These n
equations
provide the smallest error and therefore converge relatively rapidly when
compared to the
other (m-n) more slowly converging equations.
Block 135 involves explicitly solving the n equations selected in block 130 to
determine
the corresponding parameters x,, xZ,..., xõ of x. The solution x,, xz,..., xõ
of these n
equations minimizes the LAV error or cost function in (3) since we have
deliberately set
the corresponding values of v equal to zero for the n equations selected in
block 130.
Accordingly, the cost function in (3) will always be the minimum value
obtainable for
measurement equation (I) whether it is linear or nonlinear. As described in
more detail
below, it may be desirable to iterate in the nonlinear case in order to obtain
an optimal
solution.
Example Applications for Unconstrained Optimization
Example 1
We assume that we obtain the following data values in block 102 from a system
under
consideration: (0.5,1,1), (1,1,2), (1.5,1,3), (2,1,4), (2.5,1,5), (1.5,1,15),
(1.5,1,16). We
assume that the first value in each ordered triplet is a measured variable
(e.g. an output 16
which forms output vector z) and the second and third values in each ordered
triplets are
control variables (e.g. inputs 14 which form input matrix H). We assume
further that we
are trying to fit these data values to a polynomial (in this case a straight
line) of the form
CA 02620743 2008-01-10
-20-
y = a, x + a2. The polynomial y = a, x + a2 may represent the form of model
referred to
in block 56 of method 50 (Figure 2)
In block 104, we construct the following measurement equation
0.5 1 1 v,
1.0 1 2 vZ
1.5 1 3 v3
z- Hx = v= 2.0 - 1 4 az = V 4 (15)
2.5 1 5 a' v5
1.5 1 15 v6
1.5 1 16 v7
where x=[aZ,a,]". It is noted from (15) that m=7>n=2.
Block 106 involves an inquiry into whether any of the individual equations in
(15) are
identical. It can be observed that there are no identical equations in (15).
Consequently,
the block 106 inquiry is negative and method 100 proceeds to block 109. Block
109
involves normalizing (15) by dividing each of the individual equation in (15)
by a
suitable constant such that z=[z,, zZ,..., z,, ]r [0.1,0.1,...,0.1]' The block
109
normalization yields
1 1/5 1/5 v,/5
.1 1/10 2/10 v2/10
.1 1/15 3/15 v3 /15
1- 1/20 4/20 1aZ = v, /20 (15')
1 1/25 5/25 vs/25
1 1/15 15/15 v6/15
1 1/15 16/15 v7/15
Method 100 then proceeds to block 110, where we determine the a values of the
normalized matrix H in the normalized measurement equation (15'). As discussed
above,
one technique for determining these a values is method 110 of Figure 4. In the
interest of
CA 02620743 2008-01-10
-21-
brevity, we do not explicitly determine all of the a values here. However, it
may be easily
observed that
a2167 =15/15+16/15>1.
Since a 216,7 > 1(i.e. there are one or more block 110 a values greater than
or equal to
unity), the block 115 inquiry is negative and method 100 proceeds to block
117. In block
117, we divide (15') by 10 to obtain
.01 1/50 1/50 v,/50
.01 1/100 2/100 v2/100
.01 1/150 3/150 v3/150
.01 - 1/200 4/200 a2 = v4 / 200 (16)
.01 1/250 5/250 a' v5 / 250
.01 1/150 15 / 150 v6 / 150
.01 1/150 16/150 v, /150
Method 100 then proceeds to block 110 again, where we determine the a values
for the
revised equation (16). In this example, n=2 (i.e. there are two columns in the
normalized
H matrix. Consequently, there will be a ai value and an a2 value for each
combination of
2 rows. The a values calculated for the revised equation (16) are reproduced
below.
a,l1.z =laõI+la211=1/50+1/100=.03; azl,.Z =Ia12l+laz2l=1/50+2/100=.04;
a,l,,; =laõI+Ia31I=I/50+1/150=.027; aZl,.3 =Ia12I+la,Zl=1/50+3/150=.04;
a,11,4 =laõI+la411=1/50+1/200=.025; aZl1.4 =la1zl+la421=1/50+4/200=.04;
a, 1,5 =laõI+las,I=1/50+I/250=.024; a2l,,s =Ia,ZI+la521=1/50+5/250=.04;
a1Il,6 -Ia, , I +la611=1/50+1/150=.027; azl1,6 =la1zl+la621=1/50+15/150=.12;
a,11,7 = laõ I+ 1aõ I=1 / 50 + 1/ 150 =.027; a21,., =1a121 + 1a721= 1/ 50 + 16
/ 150 =.127;
a,I23=Ia21l+laõI=1/100+1/150=.017; a212.3=ja221 +la32 =2/100+3/150=.04;
a,lz,4 =1a211+1a411=1/100+1/200=.015; a2 124 =ja221 +Ia411 =2/100+4/200=.04;
a,12.5 =1a211+1a511=1/100+1/250=.014; a212.1 =ja22 +Ia521=2/100+5/250=.04;
CA 02620743 2008-01-10
-22-
a, I2,6 = laZ, l+ 1a61 l=1 / 100 + 1/ 150 =.017; aZ 12,6 = ja22 1+ 1a62 l= 2/
100 + 15 / 150 =.12;
a,I2,1 -Ia21l+la,,l=1/100+1/150=.017; azlz_7 =1a221+1a7,1=2/100+16/150 = .127;
a,13.4 =1a311+1a411=1/150+1/200=.012; az13,4 -Ia321+1a421=3/150+4/200=.04;
a,I3.s 511 =1/150+1/250=.011; aZj3,s =la3z1 +Ia52I=3/150+5/250=.04;
a113,6 -1a311 +1a611 1/150+1/150=.013; a2136 =1a321 +la621=3/150+15/150=.12;
a,13,7 =1a311+1a711=1/150+1/150=.013; a213,7 =la321+1a721=3/150+16/150=.127;
a, 14,5 =1a411 + as, j =1 / 200 + 1/ 250 =.009; a2 14,5 = 1a421 + 1a521 = 4/
200 + 5/ 250 =.04;
a,14=6 -1a411 +Ia611=1/200+1/150=.012; a214.6 =1a421+1a62,=4/200+15/150=.12;
a,14,7 =1a411+1aõ1 =1/200+1/150=.012; a214,7 =1a421 +1a721=4/200+16/150=.127;
a,I5,6 =1a51I+1a61I=1/250+1/150=.011; aZ15.6 =1a121 +la62I=5/250+15/150=.12;
a115,7 = las, I+ 1a71( = 1/ 250 + 1/ 150 = .011; a215,7 = Ia521 + 1a72 l= 5/
250 + 16 / 150 =.127;
-5
a116,7 =1a611+1a711=1/150+1/150=.013; a,16,7 =1a6,1 +la721
=15/150+16/150=.207;
It can be seen that all of the a values for the revised H matrix are strictly
less than unity.
Accordingly, the block 115 inquiry is positive and method 100 proceeds to
block 130.
Block 130 involves selecting the n equations corresponding to the row
combination with
the smallest a sum a; = a, + a z+...a ~). In this particular Example, n=2.
;-~
Consequently, block 130 involves selecting the n=2 individual equations
corresponding
to the row combination with the lowest a sum, where the a sums for a
combination of
rows j, k is detennined according to a(sum)I J k = a, I jk + aZ I j k. The a
sums
corresponding to equation (16) are given by:
a(sum)I I Z= a, 11 2+ a211 2 =.07 30 a(sum)I, 5=.064
a(sum)I, 3 = .067 a(sum)I 1 6 =.147
a(sum)I 1 4 =.065 a(sum)I I.7 =.154
CA 02620743 2008-01-10
-23-
a(sum)I Z 3 =.057 a(sum)I3 7 =.14
a(sum)I Z 4=.055 10 a(sum)I4 5=.049
a(sum)I2 5 =.054 a(sum)I4 6 =.132
a(sum)I26 =.137 a(sum)l47 =.139
a(sum)I2 7 =.144 a(sum)I5,6 =.131
a(sum)I; 4 =.052 a(sum)I S 7 =.138
a(sum)Iz,5 =.051 15 a(sum)I63 =.22
a(sum)I;6 =.133
It may be seen from the a sums calculated above that, in the case of system
(16), the
smallest a sum is a(sum)I4,5 = a114,5 + a214,s ==049 . Accordingly, in this
particular
example where n=2, block 130 involves selecting the 4`h and 5th individual
equations
20 from system (16). As discussed above, this selection of individual
equations minimizes
the cost function in (3).
Block 135 involves assuming that v4=vS=0 and explicitly solving the 4`h and
5`h
individual equations of (16)
25 O 1- 1200 ~00 az = [01 (17)
.01 1250 Y50 a' 0
Solving (17) yields aZ= (18)
a, 00.5
to provide a fitted curve of y = 0.5x.
Here it is instructive to make a number of observations about the application
of method
30 100 to this Example. Using (3), (15) and (18), we may compute the cost
function J in (3)
to be J=12.5. It can be observed that method 100 completely rejects the
outlying data
values corresponding to the 6 'h and 7'h individual equations in (15). This
can be shown by
CA 02620743 2008-01-10
-24-
letting z6=7.5 and z7=8.0 (which are the "correct" values which fit a curve of
the form
y = 0.5x ) and reapplying method 100. In such a reapplication of method 100,
the row
combination with the smallest a sum is a(sum)I6 7 = a116,7 + a 2 67 .
Accordingly, in such
a circumstance block 130 involves selecting the 6`h and 7`h individual
equations and block
135 involves solving
.01 _ Y50 Y0 a2 = 101 (19)
O1 1 g00 Y0 a` 0
aZ
which again yields = (20)
a, 00.5
This is the same solution that was obtained with the "erroneous" (e.g.
outlying) data
i0 present. However, when we calculate the cost function (3) with the new data
values, we
observe that J=O. It follows that significant errors in the measurements can
be detected
by computing the cost function J in (3). It is also possible to detect errors
in the model
curve to be fitted.
Example 2
This Example is an economic problem taken from R. Ramanathan, Introductory
Econometrics with Applications, Second Edition, Harcourt Brace Javanovich
College
Publishers, New York, 1992, pp. 135-136, which is hereby incorporated herein
by
reference and which is referred to hereinafter at "Ramanathan". This example
refers to 57
sets of data for the cumulative cost of maintenance (excluding gasoline) and
the number
of miles driven (in thousands) for a vehicle. These 57 sets of data are
obtained in block
102. The objective of method 100 in this Example is to fit the data to the
following model
cost = a,x + a2 where x= miles (in thousands).
Block 104 involves using the 57 data values to generate a measurement
equation. In this
particular example, m=57 and n=2. Blocks 106 and 108 involve removing
identical
equations, if necessary.
CA 02620743 2008-01-10
-25-
Block 109 involves normalizing the 57 individual equations. In the interest of
brevity,
this calculation is not performed explicitly here. Block 110 involves
determining the a
values for the 57 individual equations. In one particular embodiment, this may
be done in
accordance with method 110 of Figure 4. Again, in the interest of brevity,
this calculation
is not performed explicitly here. We then proceed through the loop of blocks
110, 115
and 117 until we have ensured that all of the a values are strictly less than
1.
As discussed above, block 130 involves selecting n individual equations
corresponding to
the row combination with the smallest a sum. In this Example, n=2 and the row
combination with the smallest a sum is a(sum)I44 47 = a1144,47 + a2144,41
where it may be
observed that
a1144.47 -1a44.i +1a47.11=.00395
aZ I44,47 = Ia44,2 + ja47,21 - '02g 1
Accordingly, block 130 involves selecting the 44`h and 47th individual
equations.
Block 135 involves explicitly solving the system of n=2 selected equations
(i.e. the 44`h
and 47`h individual equations). It may be seen by reference to Ramanathan that
the
solution to the system of the 44`h and 47`h individual equations is a2=
301.909
a, 55.7575
Accordingly, the fitted curve is given by CostLA1, = 55.7575x - 301.909.
Once again, it is useful to make some observations relating to this Example.
For the
solution (22), the cost function (3) may be computed to be J=34,649.6. When a
least
squares estimate (LSE) is used to evaluate this problem, on this Example, the
solution
obtained is Cost,,,,; = 53.788x - 806.08. It may be observed that the slopes
of the curves
fitted using method 100 and using a conventional LSE process are nearly equal.
When a conventional linear programming (LP) technique is used to generate a
LAV
estimate, the solution obtained is Cost,1,,.,,,,. = 18.8x -138.68 . The cost
function (3)
CA 02620743 2008-01-10
-26-
evaluated with this LAV/LP solution yields J = 44, 661.1. It may be observed
that the
fitted curves obtained from the two LAV methods do not agree.
Example 3
Consider the normalized system
0.1 0.1 0.1
z= Hx + v= 0.1 = 0.1 0.1 x' + v (17)
0.1 5 0.1 xZ
By observation it can be seen that the first two individual equations in (17)
are identical
to one another. Accordingly, when applying method 100 to equation (17), the
block 106
inquiry would be positive and method 100 would proceed to block 108. In block
108,
only one of the individual equations I and 2 would be retained. If one of
these equations
is removed, then equation (17) becomes
10.1 _ 0.1 0.1 x
z= Hx +_v = O.I 0.5 0.1 x' + v (17')
z
which actually is a system of 2 equations in 2 unknowns and therefore has a
unique
solution (i.e. it is not necessary (although it is possible) to continue with
method 100).
Example 4
Consider the following normalized system
0.1 0.1 0.1
0.1 = 0.2 0.3 x' +v (18)
0 5 0.1 IXZ
It is observed that z3=0. This represents an interesting circumstance, because
the third
individual equation can be made arbitrarily small by dividing it by any
suitable constant.
In this case, the row combination with the smallest a sum will always include
the third
row. Accordingly, the 3a individual equation in (18) will always be one of the
n=2
equations selected in block 130.
Particular Embodiment for Optimization Under Linear Constraints
Figure 5 illustrates a method 300 for optimizing the measurement equation (1)
CA 02620743 2008-01-10
-27-
z = Hx+v (1)
subject to a constraint of the form
Cx = d (19)
where C is a lxn matrix (I<n) and d is a!xl vector. It will be appreciated
that constraint
(19) actually comprises I individual equations. Constraint (19) may be the
constraint
obtained in block 58 of method 50 (Figure 2). By way of non-limiting example,
method
300 may be performed as a part of block 60 of method 50 (Figure 2).
In many respects, method 300 is similar to method 100 and blocks of method 300
corresponding to similar blocks of method 100 are preceded by the numeral "3 '
rather
than "1". Method 300 differs from method 100 in that constraint (19) must be
satisfied
exactly. In general, constraint (19) is an order-I constraint. Accordingly, we
may only
select n-I individual equations from the measurement equation (1) to be
interpolated by
the best LAV approximation Hx~~,,; This aspect of method 300 is explained in
more
detail below.
Blocks 302, 304, 306, 308 and 309 of method 300 are substantially similar to
blocks 102,
104, 106, 108 and 109 of method 100. The block 310 calculation of a values
differs from
the calculation of a values in block 110. In block I 10, each a value is
determined by a
sum of the elements of the normalized H matrix from a particular column and
from a row
combination of n rows. In block 310, however, each a value is determined by a
sum of
the elements of the normalized H matrix from a particular column and from a
row
combination of n-1 rows. This reduced number of row combinations is used to
accommodate the constraints (19).
Blocks 315 and 317 are substantially similar to corresponding blocks 115, 117
of method
100. Block 330 of method 300 differs from block 130 of method 100 in that
block 330
involves selection of the n-I individual equations from within measurement
equation (1)
which correspond to the row combination with the smallest a sum. The block 330
a sums
may be determined using substantially the same procedure as discussed above
for block
CA 02620743 2008-01-10
-28-
130. However, rather than selecting n individual equations corresponding to
the row
combination with the smallest a sum as was done in 130, block 330 involves
selecting n-1
individual equations corresponding to the row combination with the smallest a
sum.
Block 335 involves forming a system of n equations in n unknowns using the n-1
individual equations selected in block 330 together with the 1 equations from
the
constraint (19). As with block 135 of method 100, we assume that the n-1
individual
equations selected in block 330 are interpolated by the LAV approximation
Hx,.,,. . The
explicit solution to the block 335 system of n equations in n unknowns
provides the
i o individual elements of x in (1). The block 335 solution minimizes the LAV
cost function
(3) subject to constraints (19), since we have deliberately set the
corresponding values of
v equal to zero for the n-I equations selected in block 330 and for the
constraints (19).
Example Application of Optimization Under Linear Constraints
Example 5
Here we revisit the data values of Example I and we are trying to fit the
curve
y = a, x + az . However, we now make out solution subject to the constraint
that
y(10) = 4.5. We may express this constraint in the form of equation (19)
Cx = d=[1 10] a2 = 4.5 (20)
- [a,]
It is observed that 1=1 for constraint (20).
Method 300 proceeds through blocks 302, 304, 306 and 309 in manner similar to
blocks
102, 104, 106, 109 discussed above for Example 1. When method 300 arrives in
block
310 for the first time, the a values are calculated on the basis of equation
(15'). However,
in block 310, the a values are calculated using the elements from a particular
column and
from a row combination of n-1=2-1=1 row. For the sake of brevity, we do not
perform
this calculation here. However, it is observed that
=1
a116 = Ja621
and
CA 02620743 2008-01-10
-29-
a117 =Ia721=16/15
Accordingly, the block 315 inquiry is negative and method 300 loops through
block 317
where the measurement equation (15') is divided by 10 to obtain the
measurement
equation (16). Method 300 then retums to block 310, where the a values are
calculated
again. The a values calculated in block 310 are:
a,11 =laõI=1/50=.02; a211 =Ia121=1/50=.02;
a112 =1a211=1/100=.01; a212 =1a221 =2/100=.02;
a1I3 -Ia;fl=1/150=.007; a213 -la?21-3/150=.02;
a,14 =1a411 =1/200=.005; a214 =1a421=4/200=.02;
a, Is ='as, I=1 / 250 =.004; a2'5 = ja52 1= 5/ 250 = .02;
a116 -1a611 1/150=.007; a216 =ja621=15/150=.1;
a1 17 =1a71 l =1 / 150 = .007; a217 = la,21 =16 / 150 = .107;
Since all the a values are less than unity, the block 315 inquiry is positive
and method
300 proceeds to block 330. Block 330 involves selecting the n-1=2-1=1 equation
corresponding to the row combination with the lowest a sum. In this Example,
the row
"combinations' are only single rows and the single row with the lowest a sum
is the 5`n
row which can be observed to have an a sum of a(sum)I s= a, Is + a215 =.024.
Accordingly, in this Example, the individual equation of the 5`h row is
selected in block
330.
In block 335, we form a system of equations using the constraint (20) and the
individual
equation(s) selected in block 330 (i.e. the equation corresponding to the 5`h
row of (16))
under the assumption that vi in (1) is zero for the individual equation(s)
selected in block
130. For this Example, the resultant system is given by:
1 10 a2 4.5 (21)
1/250 5/250 al .01
CA 02620743 2008-01-10
-30-
which may be solved to yield a2 =.4 . This solution, which corresponds to
i
y =.4x +.5, differs from the solution calculated in Example 1. However, this
solution
minimizes the cost function (3) subject to the constraint (20).
Non-Linear Optimization
In other embodiments, methods similar to those discussed above may be applied
to non-
linear measurement equations of the form
z = f (a, x) + v (22)
where z represents an mxl measurement vector (e.g. from output 16 (Figure 1)),
x is an
nxl vector representing the parameters x,, xZ,..., xõ of the curve/model (e.g.
model 26 of
system 12 (Figure 1)) where m>n, a represents the explanatory variables (e.g.
from inputs
14 (Figure 1)) and v represents the measurement error. Measurement equation
(22) may
be expressed in component form as
z; = f,.(a;,x)+v, i= l, 2===m (23)
In the non-linear case (22), the LAV cost function is given by
m m
J = E I zf -.f (ai,x)I = I Iv, l (24)
;=1 ;_1
We also impose a constraint (expressed in component form)
g,(ar,x) = 0 (25)
Constraint (25) may be the constraint obtained in block 58 of method 50
(Figure 2), for
example..
We may linearize the component form measurement equation (23) by first order
Taylor
series expansion about an initial estimate xo to provide
8f, (a;, x)
z; = f,~(aõx,)+v, + (x-x ) (26)
C7X s=z
We can then define the following
Ox = (x - x,F) (27)
CA 02620743 2008-01-10
-31-
Az; =z; - f,(a;,& ) fori=1,2, ... m (28)
H; = (a~'X) for i=1,2, ... m (29)
s-~
and rewrite equation (26) as
Az, = H;Ox + v; for i=1,2, ... m (30)
or, in complete matrix/vector form,
Oz = HOx + v (31)
In a similar manner, the constraint (25) may be linearized by Taylor series
expansion to
provide
C,Ax = d, for i=1,2, ... 1 (32)
or, in complete matrix/vector form,
CAx = d (33)
where
C= ag'(a''x) for i=1,2, ... 1 (34)
ax x-x
d, _ -g;(a,,xõ) for i=1,2, ... 1 (35)
Examining (31) and (33)
Az = HOx + v (31)
COx = d (33)
it may be recognized that these equations have the same form as the
optimization under
linear constraints presented above in equations (1) and (19), except that z
and x of
equations (1) and (19) are replaced by dz and dx in equations (31) and (33).
Equations (31) and (33) may therefore be optimized according to the procedure
presented
above for optimization under linear constraints to yield a solution dx. The
original
estimated solution xo may then be iterated by increments of dx of to solve the
original
non-linear problem posed by equations (23) and (25). More particularly, the
(i+1)'ti
iterate is given by
CA 02620743 2008-01-10
-32-
x;,, = x; + Ox, where i=iteration counter (36)
The iteration continues until dx; is sufficiently small to meet applicable
acceptance
criteria. For example, the iteration may continue until
Ax; < 0 (37)
where 0 represents an iteration threshold. In some embodiment, iteration
criteria may
comprise other thresholds. By way of non-limiting example, iteration criteria
may be
based on threshold value of the cost function (24), a number of iterations or
the like.
Figure 6 represents a method 400 for optimizing a non-linear system according
to a
t o particular embodiment of the invention. By way of non-limiting example,
method 400
may be implemented in block 60 of method 50 (Figure 2). Method 400 commences
in
block 402. Block 402 involves obtaining data values (e.g. from inputs 14 and
outputs 16)
and is similar to block 102 of method 100 discussed above. Block 404 involves
using the
block 402 data values to generate a non-linear measurement equation of the
form of (23)
and a constraint equation of the form of (25).
Method 400 then proceeds to block 410 which involves linearizing the block 404
measurement and constraint equations in the manner described above to arrive
at a linear
optimization problem in dx - i.e. in the form described above in equations
(31) and (33).
Method 400 then proceeds to block 415 which involves making an initial
estimate of the
value of x in the block 404 measurement equation. In some applications, there
may be
some basis on which to form the estimate of the parameters xl, X2, ... x, such
as past
experiment, knowledge of system 12 under consideration or the like. Block 415
typically
involves substituting the initial estimate into the block 410 linearized
measurement and
constraint equations to arrive at estimated linearized measurement and
constraint
equations.
Method 400 then proceeds to block 420. Block 420 involves carrying out the
steps of
blocks 306 to 335 of method 300 on the estimated linearized measurement and
constraint
3o equations from block 415. The application of blocks 306 to 335 may be
substantially
similar to the application discussed above and the result of block 420 is to
determine a
CA 02620743 2008-01-10
-33-
value for Ax based on the estimated linearized measurement and constraint
equations
from block 415. In block 425, the block 420 Ax value is used to update the
iterated
estimate of x according to equation (36).
Block 430 involves an inquiry as to whether the block 420 Ax satisfies
suitably selected
convergence criteria. If the block 420 dx satisfies the convergence criteria
(block 430
YES output), then method 400 proceeds to block 435 where the block 425
estimate ofx is
provided as the final estimate of x- e.g. the parameters x,, x2 ,..., xn of
model 26 (Figure
1). If the block 420 Ax does not satisfy the convergence criteria (block 430
NO output),
lo then method 400 proceeds to block 440 where the linearized measurement and
conversion equations are updated with the new estimate of x (i.e. as updated
in block
425) before looping back to block 420, where a new Ax value is obtained.
Method 400
loops through blocks 420, 425, 430 and 440 until the block 430 convergence
criteria are
satisfied. Typically the block 430 convergence criteria include threshold
conditions
relating to one or more parameters, such as Ax, the LAV cost function (24),
the number of
iterations or the like must be less than the threshold.
Figure 7 schematically depicts a method 500 for optimizing a non-linear system
according to another embodiment of the invention. By way of non-limiting
example,
method 500 may be implemented in block 60 of method 50 (Figure 2). Method 500
is
similar in many respects of method 400 and blocks of method 500 corresponding
to
similar blocks of method 400 are preceded by the numeral "5" rather than "4".
Blocks
502, 504, 510, 515 and 520 of method 500 are substantially similar to
corresponding
blocks 402, 404, 410, 415 and 420 of method 400.
Method 500 differs from method 400 in that, after the first iteration, method
500 proceeds
on the assumption that the n-1 individual equations selected in block 520
(more
particularly, in the block 330 aspect of block 520) will be the n-1 selected
equations
throughout the iteration. Method 500 can therefore save computational
resources by only
recalculating necessary components of the measurement equation and convergence
criteria for each iteration. This aspect of method 500 is explained in more
detail below.
CA 02620743 2008-01-10
-34-
After determining a first a new dx in block 520, method 500 proceeds to the
block 530
inquiry as to whether the block 520 dx values satisfy suitably selected
convergence
criteria. The block 530 inquiry is substantially similar to the block 430
inquiry described
above. If the block 530 convergence criteria are not satisfied (block 530 NO
output), then
method 500 proceeds to block 525. Block 525 is substantially similar to block
425 of
method 400 and involves updating the iterated estimate of x according to
equation (36)
using the Ax obtained in the previous iteration.
Method 500 then proceeds to block 540. Block 540 involves using the updated
iterated
estimate of x from block 525 to update the n-1 individual linearized
measurement
equations selected in block 520 (more particularly, in the block 330 aspect of
block 520).
Method 500 assumes that the n-1 individual linearized measurement equations
selected in
block 520 (i.e. in the first iteration) will be the same in each iteration.
Consequently, it is
not necessary to update the m-(n-1) other, non-selected individual equations
in the
linearized measurement equation for the second and subsequent iterations. This
aspect of
method 500 conserves computational resources. Block 540 also involves using
the
updated iterated estimate of x from block 525 to update the 1 individual
linearized
constraint equations.
The result of updating the n-I individual linearized measurement equations and
the 1
individual linearized constraint equations with the new updated iterated
estimate of x
from block 425 is a system of n equations in n unknowns dxl, dx2, ... Ax,,.
Block 540 also
involves solving this system of equations to arrive at a new value for Ax.
Method 500
then returns to the block 530 inquiry to determine whether the block 540
values of Ax
satisfy the convergence criteria. Method 500 continues to loop through blocks
530, 525
and 540 until the convergence criteria are satisfied in block 530.
When the block 530 convergence criteria are satisfied (block 530 YES output),
method
3o 500 proceeds to block 532. Block 532 is similar to block 525 and involves
updating the
iterated estimate of x according to equation (36) using the dx obtained in the
previous
CA 02620743 2008-01-10
-35-
iteration. In block 535, the block 532 value of x is presented as the final
estimated
solution - e.g. the parameters x,, x2,..., xõ of mode126 (Figure 1).
Example Application of Non-Linear Optimization
Example 6
We consider the non-linear relations
f,(x)= f,(x,,x,)=x; +xz-10
f2(x)=f2(x11x2)=x1 +x2 -7 (38)
and we desire to minimize F(x) = I f(xA + I fZ (X_)I (39)
subject to the constraint g, (x, , x2 )= x; - x2 -1= 0 (40)
In block 404, we use (38) to generate a measurement equation of the form of
equation
(23)
z, = f,(x)+v, = f,(x, ,x,)+v, = x; +xZ -10+v1
z2 = fZ(x)+v2 = f,(x,,xz)+vz =x, +xZ -7+v2 (41)
Comparing (39) to the general cost function (24) described above, it may be
seen that
minimizing (39) is equivalent to minimizing the general cost function (24)
when z;=0.
Accordingly, in block 410, we set z;=0 in the measurement equation (41) and
then
linearize the resultant equations according to the procedure described above
in equations
(27), (28), (29) and (31). This linearization process yields
Oz = HAx + v (42)
where
Az - (x; + xz -10)
- (x, + x2 - 7) (43)
H = 2x , 1 (44)
1 2x2
Ax = ~' (45)
z
CA 02620743 2008-01-10
-36-
Block 410 also involves linearizing the constraint (40) using the process
described above
in equations (33), (34) and (35) to obtain
CAx = d (46)
where
C =12x, - 3x2 j (47)
d = [-(x; - x2 -1)] (48)
Block 415 involves making an initial estimate of xo. In this Example, we
assume
x
xo = x'" = 1 and substituting these xo values into the linearized equations
(42)-(48).
zõ
Upon substituting the xo estimate, equation (42) becomes
Oz=HAx+v= 8 - 2 1 Axj +v (49)
- - - - 5 1 2 Ax2 -
and equation (46) becomes
CAx=d=[2 -3 ~' =1 (50)
z
Method 400 then proceeds to block 420 which involves using the procedures
outlined in
blocks 306 to 335 of method 300 to solve for Ax = Ox~' . For the sake of
brevity, we do
_ z
not describe all of the individual blocks of method 300. However, it may be
seen that,
once normalized, equation (49) becomes
.l .025 .0125 dz,
Oz=HOx+v= 1 .02 .04 ~ +v (51)
- - z
The a values for (51) include row combinations of n-1 rows. In this case, n-
1=2-1=1 (i.e.
row "combinations" of a single row). The a values and may be calculated as
a,l, =aõI=.025; aZl1 =a12 =.0125;
a 112 =laz,I=.02; a 212 =lazzl=.04;
CA 02620743 2008-01-10
-37-
The lowest a sum is a(sum)I , = a, I, + a z I, =.0375. We therefore select the
first
individual equation from (51) together with the constraint (50) to develop a
system of
n=2 equations in n=2 unknowns which is given by
025 .0125 Ox, .1
2 -3 OxZ 1 (52)
which yields [AXI 3.125 AX= 1 75 (53)
z
Method 400 then proceeds to block 425 which involves updating the iterated
estimate
according to equation (36). In this case
x, _ x, + dx, _ 1+ 3.125 _ 4.125 (54)
xz xz o OzZ 1 1.75 2.75
In block 430, method considers whether or not the convergence criteria are
met. In this
particular Example, a non-limiting example of a suitable convergence criteria
may be
dxI, dz2<0.05. Since this convergence criteria is not satisfied (block 430 NO
output),
method 400 proceeds to block 440. Block 440 involves updating the linearized
equations
(42)-(48) with the new x estimate from (54).
Method 400 then loops back to block 420 for a subsequent iteration. The
iterations of
blocks 420, 425, 430, 440 are repeated until the convergence criteria are met
in block
430. For the sake of brevity, the 2 a and subsequent iterations are not
explicitly calculated
here. Others have shown (using different techniques) that this problem
converges to
approximately x, = 2.8425
x2 1.9202
Other Considerations
Certain implementations of the invention comprise computer processors which
execute
software instructions which cause the processors to perform a method of the
invention.
The invention may also be provided in the form of a program product. The
program
CA 02620743 2008-01-10
-38-
product may comprise any medium which carries a set of computer-readable
instructions
which, when executed by a data processor, cause the data processor to execute
a method
of the invention. The program product may be in any of a wide variety of
forms. The
program product may comprise, for example, physical media such as magnetic
data
storage media including floppy diskettes, hard disk drives, optical data
storage media
including CD ROMs, DVDs, electronic data storage media including ROMs, flash
RAM,
or the like. The program product may also comprise data, databases or other
information
which may be accessible to, but not necessarily executed by, a processor.
io Where a component (e.g. a software module, processor, assembly, device,
circuit, etc.) is
referred to above, unless otherwise indicated, reference to that component
(including a
reference to a "means") should be interpreted as including as equivalents of
that
component any component which performs the function of the described component
(i.e.,
that is functionally equivalent), including components which are not
structurally
equivalent to the disclosed structure which performs the function in the
illustrated
exemplary embodiments of the invention.
As will be apparent to those skilled in the art in the light of the foregoing
disclosure,
many alterations and modifications are possible in the practice of this
invention without
departing from the spirit or scope thereof. For example:
= The methods described above and shown in the associated Figures represent
particular embodiments of the invention. It will be appreciated that the
procedures
associated with various blocks may be performed in different orders than the
one
shown herein. By way of non-limiting example, in method 50 (Figure 2), the
model to be considered (block 56) and/or the constraints to be considered
(block
58) may be obtained prior to obtaining inputs and outputs (blocks 52, 54). In
addition, the temporal performance of the procedures associated with various
blocks may overlap. By way of non-limiting example, in method 100 (Figure 3)
and similar methods, it may be possible to check if a values are strictly less
than
unity (block 115) at a time that overlaps with determination of the a values
(block
110). Such a temporal overlap may be desirable because if it is discovered
that
CA 02620743 2008-01-10
-39-
any one a value is greater than or equal to unity, then it may not be
necessary to
compute further a values.
= Once a model characterized by a nxl parameter vector x=[xI, X2, ... x,,JT is
obtained using one of the methods described herein, the model may be used to
predict outputs 16 of system 12 under consideration based on a known set of
inputs 14.
As such, the invention should be interpreted in accordance with the claims
appended
hereto.
Appendix A
In [14] pp 45-46 the authors pose the following question: Under what
conditions is the
mapping A a contraction when it is defined by Ax = x? Then the following
condition is
derived.
Let d(x, y) =ZI x; - y; l then
i_1
d(Y',Y") - Y '-Y ~~ ) a (x x ") S a . x .'-x
-~iI ~ + I-~i L~j 9 .l I ~;~jl iJ I J l
5 max j Y
,, l a;j I d(x', x" ) Hence we have the following condition for contraction:
Zlayl<_a<I.
