1

A New Continuous-Time
Equality-Constrained Optimization Method to
Avoid Singularity

arXiv:1209.5218v2 [] 8 Jan 2014

Quan Quan and Kai-Yuan Cai

Abstract
In equality-constrained optimization, a standard regularity assumption is often associated with
feasible point methods, namely the gradients of constraints are linearly independent. In practice, the
regularity assumption may be violated. To avoid such a singularity, we propose a new projection
matrix, based on which a feasible point method for the continuous-time, equality-constrained optimization problem is developed. First, the equality constraint is transformed into a continuous-time
dynamical system with solutions that always satisfy the equality constraint. Then, the singularity
is explained in detail and a new projection matrix is proposed to avoid singularity. An update (or
say a controller) is subsequently designed to decrease the objective function along the solutions of
the transformed system. The invariance principle is applied to analyze the behavior of the solution.
We also propose a modified approach for addressing cases in which solutions do not satisfy the
equality constraint. Finally, the proposed optimization approaches are applied to two examples to
demonstrate its effectiveness.

Index Terms
Optimization, equality constraints, continuous-time dynamical systems, singularity

I. I NTRODUCTION
According to the implementation of a differential equation, most approaches to continuoustime optimization can be classified as either a dynamical system [1],[2],[3] or a neural
network [4],[5],[6],[7]. The dynamical system approach relies on the numerical integration
of differential equations on a digital computer. Unlike discrete optimazation methods, the
The authors are with Department of Automatic Control, Beijing University of Aeronautics and Astronautics, Beijing
100191, China (Quan Quan: qq buaa@buaa.edu.cn, http://quanquan.buaa.edu.cn/; Kai-Yuan Cai: kycai@buaa.edu.cn)
January 9, 2014

DRAFT

2

step sizes of dynamical system approaches can be controlled automatically in the integration
process and can sometimes be made larger than usual. This advantage suggests that the
dynamical system approach can in fact be comparable with currently available conventional
discrete optimal methods and facilitate faster convergence [1],[3]. The application of a higherorder numerical integration process also enables us to avoid the zigzagging phenomenon,
which is often encountered in typical linear extrapolation methods [1]. On the other hand,
the neural network approach emphasizes implementation by analog circuits, very large scale
integration, and optical technologies [8]. The major breakthrough of this approach is attributed
to the seminal work of Hopfield, who introduced an artificial neural network to solve the
traveling salesman problem (TSP) [9]. By employing analog hardware, the neural network
approach offers low computational complexity and is suitable for parallel implementation.
For continuous-time equality-constrained optimization, existing methods can be classified
into three categories [1]: feasible point method (or primal method), augmented function
method (or penalty function method), and the Lagrangian multiplier method. Determining
whether one method outperforms the others is difficult because each method possesses distinct
advantages and disadvantages. Readers can refer to [1],[4],[7],[10] and the references therein
for details. The feasible point method directly solves the original problem by searching
through the feasible region for the optimal solution. Each point in the process is feasible,
and the value of the objective function constantly decreases. Compared with the two other
methods, the feasible point method offers three significant advantages that highlight its
usefulness as a general procedure that is applicable to almost all nonlinear programming
problems [10, p. 360]: i) the terminating point is feasible if the process is terminated before
the solution is reached; ii) the limit point of the convergent sequence of solutions must
be at least a local constrained minimum; and iii) the approach is applicable to general
nonlinear programming problems because it does not rely on special problem structures
such as convexity.
In this paper, a continuous-time feasible point approach is proposed for equality-constrained
optimization. First, the equality constraint is transformed into a continuous-time dynamical
system with solutions that always satisfy the equality constraint. Then, the singularity is
explained in detail and a new projection matrix is proposed to avoid singularity. An update
(or say a controller) is subsequently designed to decrease the objective function along the
solutions of the transformed system. The invariance principle is applied to analyze the
behavior of the solution. We also propose a modified approach for addressing cases in which
January 9, 2014

DRAFT

3

solutions do not satisfy the equality constraint. Finally, the proposed optimization approach
is applied to two examples to demonstrate its effectiveness.
Local convergence results do not assume convexity in the optimization problem to be
solved. Compared with global optimization methods, local optimization methods are still
necessary. First, they often server as a basic component for some global optimizations, such
as the branch and bound method [11]. On the other hand, they can require less computation
for online optimization. Compared with the discrete optimal methods offered by MATLAB,
at least two illustrative examples show that the proposed approach avoids convergence to a
singular point and facilitates faster convergence through numerical integration on a digital
computer. In view of these, the contributions of this paper are clear and listed as follows.
i) A new projection matrix is proposed to remove a standard regularity assumption that
is often associated with feasible point methods, namely that the gradients of constraints are
linearly independent, see [1, p.158, Equ.(4)],[2, p.156, Equ.(2.3)],[7, p.1669, Assumption 1].
Compared with a commonly-used modified projection matrix, the proposed projection matrix
has better precision. Moreover, its recursive form can be implemented more easily.
ii) Based on the proposed matrix, a continuous-time, equality-constrained optimization
method is developed to avoid convergence to a singular point. The invariance principle is
applied to analyze the behavior of the solution.
iii) The modified version of the proposed optimization is further developed to address cases
in which solutions do not satisfy the equality constraint. This ensures its robustness against
uncertainties caused by numerical error or realization by analog hardware.
We use the following notation. Rn is Euclidean space of dimension n. k·k denotes the
Euclidean vector norm or induced matrix norm. In is the identity matrix with dimension n.
0n1 ×n2 denotes a zero vector or a zero matrix with dimension n1 × n2 . Direct product ⊗
and vec (·) operation are defined in Appendix A. The function [·]× : R3 → R3×3 with matrix
H ∈ R9×3 is defined in Appendix B. Suppose g : Rn → R. The gradient of the function
g is given by ∇g (x) = ∇x g (x) = [∂g (x) /∂x1 · · · ∂g (x) /∂xn ]T ∈ Rn and the matrix
of second partial derivatives of g (x) known as Hessian is given by ∇xx : R → Rn×n and
∇xx g (x) = [∂ 2 g (x) /∂xi ∂xj ]ij .

January 9, 2014

DRAFT

4

II. P ROBLEM F ORMULATION
A. Equality-Constrained Optimization
The class of equality-constrained optimization problems considered here is defined as
follows:
min v (x) , s.t. c (x) = 0

(1)

x∈Rn

where v : Rn → R is the objective function and c = [c1 c2 · · · cm ]T ∈ Rm , ci : Rn →
R are the equality constraints. They are both twice continuously differentiable. Denote by
h
i
∇c (x) , ∇c1 (x) ∇c2 (x) · · · ∇cm (x) ∈ Rn×m . To avoid a trivial case, suppose the

constraint (or feasible set) F = { x ∈ Rn | c (x) = 0} =
6 ∅.

Definition 1 [12, pp. 316-317]. For the problem (1), a vector x∗ ∈ F is a global minimum
if v (x∗ ) ≤ v (x) , ∀x ∈ F ; a vector x∗ ∈ F is a local (strict local) minimum if there is a
neighborhood N of x∗ such that v (x∗ ) ≤ v (x) (v (x∗ ) < v (x)) for x ∈ N ∩ F .
Definition 2 [10, p. 325]. A vector x∗ ∈ F is said to be a regular point if the gradient
vectors ∇c1 (x∗ ) , ∇c2 (x∗ ) , · · · , ∇cm (x∗ ) are linearly independent. Otherwise, it is called a
singular point.
This paper aims to propose an approach to continuous-time, equality-constrained optimization to identify the local minima based on a feedback control perspective.
Remark 1. Inequality-constrained optimizations can be transformed into equality-constrained
optimizations by introducing new variables. For example, the inequality constraint x ≤ 1, x ∈
R can be replaced with an equality constraint x+z 2 = 1, z ∈ R. Also, the inequality constraint
−1 ≤ x ≤ 1, x ∈ R can be replaced with an equality constraint x = sin (z) , z ∈ R. Here,
we only focus on equality-constrained optimization.
B. Equality Constraint Transformation
Optimization problems are often solved by using numerical iterative methods. For an
equality-constrained optimization problem, the major difficulty lies in ensuring that each
iteration satisfies the constraint and can further move toward the minimum. To address this
difficulty, a transformation of the equality constraint is proposed, which is formulated as an
assumption.
Assumption 1. For a given x0 ∈ F , there exists a function f : Rn → Rn×l such that
ẋ (t) = f (x (t)) u (t) , x (0) = x0

January 9, 2014

(2)

DRAFT

5

with solutions that satisfy x (t) ∈ Fu (x0 ) , where Fu (x0 ) = {x (t) ∈ F |ẋ (t) = f (x (t)) u (t) ,
x (0) = x0 ∈ F , ∀u (t) ∈ Rl , t ≥ 0}.
From a feedback control perspective, the update u can be considered as a control input.
The objective function v (x) can be considered a Lyapunov-like function, although v(x) is
not required to be a Lyapunov function. Based on Assumption 1, the objective of this paper
can be restated as: to design a control input u to decrease v(x) along the solutions of (2)
until x has achieved a local minimum. In the following, we will omit the variable t except
when necessary.
Remark 2. The proposition of Assumption 1 is motivated by the property of attitude
kinematics [13, p. 200]: ẋ = 12 E (x) w, where x = [q0 q T ]T ∈ R4 , q0 ∈ R, q, w ∈ R3 and
E (x) = [−q q0 I3 + [q]T× ]T ∈ R4×3 . The function [·]× : R3 → R3×3 is defined in Appendix B.
All solutions of the attitude kinematics satisfy the constraint kxk2 = 1 driven by any w ∈ R3 .
The explanation is given as follows. It is easy to check that xT ẋ = 12 xT E (x) w = 0 since
[q]× q = 0 for ∀q ∈ R3 . Therefore, the solution always satisfies the constraint kx (t)k2 = 1 if
kx (0)k = 1, t ≥ 0. Another representation of attitude kinematics is
Ṙ = [w]× R

(3)

where R ∈ R3×3 is a rotation matrix satisfying the constraint RT R = I3 . For (3), we have

d
RT R = RT Ṙ + ṘT R
dt


= RT [w]× + [w]T× R = 03×3 .

That is why the evolution of R always lies on the constraint RT R = I3 .
Remark 3. The best choice of f (x) is to satisfy Fu (x0 ) = F . However, it is difficult
to achieve. For example, if c (x) = (x1 + 1) (x1 − 1), x = [x1 x2 ]T ∈ R2 , then F =
{ x ∈ R2 | x1 = 1, x1 = −1}. Since the two sets {x ∈ R2 | x1 = 1} and {x ∈ R2 | x1 = −1} are
not connected, the solution of (2) starting from either set cannot access the other. Although
Fu (x0 ) 6= F , we still expect the global minimum x∗ ∈ Fu (x0 ) . That is why we often require
that the initial value x0 be close to the global minimum x∗ . Besides this, it is also expected
that the function f (x) is chosen to make the set Fu (x0 ) as large as possible so that the
probability of x∗ ∈ Fu (x0 ) is higher.
If c (x) = Ax, A ∈ Rm×n , then the function f (x) can be chosen to satisfy F = Fu (x0 ) ,
∀x0 ∈ F .

January 9, 2014

DRAFT

6

Theorem 1. Suppose that c (x) = Ax and f (x) = A⊥ , where A⊥ is with full column rank,
and the space spanned by the columns of A⊥ is the null space of A. Then F = Fu (x0 ) ,
∀x0 ∈ F .
Proof. Since Fu (x0 ) ⊆ F , the remaining task is to prove F ⊆ Fu (x0 ) , ∀x0 ∈ F , namely
for any x̄ ∈ F there exists a control input u ∈ Rl that can transfer any initial state x0 ∈ F
to x̄. Since x0 , x̄ ∈ F , there exist u0 , ū ∈ Rl such that x̄ = A⊥ ū and x (0) = A⊥ u0 by the
definition of A⊥ . Design a control input

 1 (ū − u ) , 0 ≤ t ≤ t̄
0
t̄
.
u (t) =

0,
t > t̄.
With the control input above, we have
Z t
x (t) − x (0) =
A⊥ u (s) ds
0

=

Z

t̄

A⊥ u (s) ds = A⊥ ū − A⊥ u0 ,

0

when t ≥ t̄. Then x (t) = x̄, t ≥ t̄. Hence F ⊆ Fu (x0 ) , ∀x0 ∈ F . Consequently, F =
Fu (x0 ) , ∀x0 ∈ F . 
From the proof of Theorem 1, the choice of f (x) becomes a controllability problem.
However, it is difficult to obtain a controllability condition of a general nonlinear system.
Correspondingly, it is difficult to choose f (x) for a general nonlinear function c (x) to satisfy
F = Fu (x0 ) . Motivated by the linear case above, we aim to design a function f (x) whose
range is the null space of ∇c (x)T for any fixed x ∈ Rn . This idea can be formulated as
V1 (x) = V2 (x), where
V1 (x) = {z ∈ Rn |∇c (x)T z = 0},
V2 (x) = {z ∈ Rn |z = f (x) u, u ∈ Rl }.
III. S INGULARITY

AND

A N EW P ROJECTION M ATRIX

A. Singularity
The function f is the projection matrix, which orthogonally projects a vector onto the null
space of ∇cT . One well-known projection matrix is given as follows [1],[2],[7]:

−1 T 
f (x) = In − ∇c ∇cT ∇c
∇c (x) .

January 9, 2014

(4)

DRAFT

7

We can easily verify that ∇c (x)T f (x) ≡ 0. This projection matrix requires that ∇c (x)
should have full column rank, i.e., every x ∈ F is a regular point. However, the assumption does not hold in cases where ∇c (x)T ∇c (x) is singular. This condition is the major
motivation of this paper. For example, consider an equality constraint as
c (x) = (x1 − x2 + 2) (x1 + x2 ) = 0,
h

iT

∈ R2 . The feasible set is either { x ∈ R2 | x1 − x2 + 2 = 0} or
h
iT
{ x ∈ R2 | x1 + x2 = 0} . As shown in Fig.1, the point xp1 = −2 0
has a unique feasible
h
iT
direction and the point xp2 = 0 0
also has a unique feasible direction. Whereas, the
h
iT
point xp3 = −1 1
has two feasible directions. This causes the singular phenomena.

where x =

x1 x2

The singularity often occurs at the intersection of the feasible sets, where exist non-unique
feasible directions. Mathematically, ∇c (x)T ∇c (x) is singular. Concretely, the gradient vector
of c (x) is



∇c (x) = 

2x1 + 2
−2x2 + 2



.

At the points xp1 and xp2 , the gradient vector of c (x) is


 
−2
2
 , ∇c (xp2 ) =  
∇c (xp1 ) = 
2
2

and by (4), the projection matrices are further




0 1
0 −1
 , f (xp2 ) = 

f (xp1 ) = 
1 0
−1 0

respectively. Whereas, at the point xp3 , the gradient vector of c (xp3 ) is
 
0
∇c (xp3 ) =   .
0
−1

For such a case, ∇c (xp3 )T ∇c (xp3 )
does not exist.

To avoid singularity, a commonly-used modified projection matrix is given as follows

−1 T 
T
f (x) = In − ∇c εIm + ∇c ∇c
∇c (x)
(5)

where ε > 0 is a small positive scale. We have ∇c (x)T f (x) 6= 0 no matter how small ε is.
On the other hand, to obtain f (x) by (5), a very small ε will cause ill-conditioning problem

January 9, 2014

DRAFT

8

x1  x2  2 0

x2

> 1 1@

T

x p3

>0

x p2
x p1

> 2

0@

T

x1

0@

T

x1  x2

Fig. 1.

0

Singularity Example

especially for a low-precision processor. For example, consider the following gradient vectors:
h
i
∇c1 = 1 1 1 1
h
i
∇c2 = 2 1 1 1
h
i
(6)
∇c3 = 3 2 2 2 .
as the precision error, we employ (5) with different ε = 10−k , k =

Taking ep = ∇cT f

1, · · · , 15 to obtain the projection matrix f . As shown in Fig.2, the error varies with different
k. The best precision error can be achieved only at ε = 10−8 with a precision error around
10−8 . Reducing ε further will increase the numerical error.
0

10

−1

10

−2

10

−3

ep

10

−4

10

−5

10

−6

10

−7

10

−8

10

0

5

10

15

k

Fig. 2.

Precision error of a common-used modified projection matrix with different ε = 10−k

January 9, 2014

DRAFT

9

The best cure is to remove the linearly dependent vector directly from ∇c (x). For example,
h
i
in ∇c (x) = ∇c1 (x) ∇c2 (x) ∇c3 (x) ∈ Rn×3 , if ∇c3 (x) can be represented by a linear
combination of ∇c1 (x) and ∇c2 (x) , then ∇c (x)T ∇c (x) is singular. The best cure is to

remove ∇c3 (x) from ∇c (x), resulting in
h
i
∇cnew (x) = ∇c1 (x) ∇c2 (x) ∈ Rn×2 .

With it, the projection matrix becomes

−1 T 
T
fnew (x) = In − ∇cnew ∇cnew ∇cnew
∇cnew (x) .

It is easy to see that ∇c (x)T fnew (x) ≡ 0. For a linear time-invariant matrix ∇c (x) , namely

independent of x, we can avoid singularity by removing dependent terms out of ∇c (x)
before computing a projection matrix. However, this idea does not work for a general ∇c (x)
depending on x. Therefore, “the best cure” cannot be implemented continuously, which
further cannot be realized by analog hardware. For such a purpose, we will propose a new
projection matrix.
B. A New Projection Matrix
For a special case c : Rn → R, such a f (x) is designed in Theorem 2. Consequently, a
method is proposed to construct a projection matrix for a general case c : Rn → Rm . Before
the design, we have the following preliminary results.
Lemma 1. Let
W1 = {z ∈ Rn |LT z = 0}
T

W2 = {z ∈ Rn |z =

 1
n
where L ∈ R and δ (x) =
 0

In −

LL

δ kLk2 + kLk2

x = 0, x ∈ R
x 6= 0, x ∈ R

!

u, u ∈ Rn },

. Then W1 = W2 .

Proof. See Appendix C. 

Theorem 2. Suppose that c : Rn → R and the function f (x) is designed to be
∇c (x) ∇c (x)T

.
f (x) = In −
δ k∇c (x)k2 + k∇c (x)k2

(7)

Then Assumption 1 is satisfied with u ∈ Rn and V1 (x) = V2 (x) .

Proof. Since ċ (x) = ∇c (x)T ẋ and ẋ = f (x) u, the function f (x) is defined as in (7) so
that ċ (x) ≡ 0 by Lemma 1. Therefore, Assumption 1 is satisfied with u ∈ Rn . Further by
Lemma 1, V1 (x) = V2 (x) . 
January 9, 2014

DRAFT

10

Theorem 3. Suppose that c : Rn → Rm and the function f (x) is in a recursive form as
follows:
f0 = In


fk = fk−1 In − 
δ



T
fk−1
∇ck ∇cTk fk−1
,

2
2
T
T
∇ck
+ fk−1
∇ck
fk−1

(8)

k = 1, · · · , m. Then Assumption 1 is satisfied with f = fm and u ∈ Rn and V1 (x) = V2 (x) .
Proof. See Appendix D. 
T
∇ck 6= 0, then δ
Remark 4. In (8), if fk−1

fk = fk−1

In −



T
fk−1
∇ck

2



T
fk−1
∇ck ∇cTk fk−1
T
∇ck
fk−1

2

= 0, namely
!
.

This is the normal way to construct a projection matrix. On the other hand, if ∇ck can
T
T
be represented by a linear combination of ∇ci , then fk−1
∇ck = 0 as fk−1
∇ci = 0, i =


2
T
1, · · · , k − 1. In this case, δ fk−1
∇ck
6 0. Consequently, the projection matrix will
=

reduce to the previous one fk = fk−1 , that is equivalent to removing the term ∇ck . This is
consistent with “the best way”.

Remark 5. In practice, the impulse function δ (x) is approximated by some continuous
functions such as δ (x) ≈ e−γ|x| , where γ is a large positive scale. Let us revisit the example
for the gradient vectors (6). Taking ep = ∇cT f

as the error again, we employ (8) with

γ = 30 to obtain the projection matrix f with ep = 2.7629 ∗ 10−10 . This demonstrates the
advantage of our proposed projection matrix over (5). Furthermore, compared with (4) or (5),
the explicit recursive form of the proposed projection matrix is also easier for the designer
to implement.
IV. U PDATE D ESIGN

AND

C ONVERGENCE A NALYSIS

In this section, by using Lyapunov’s method, the update (or say controller) u is designed
to result in v̇ (x) ≤ 0. However, the objective function v (x) is not required to be positive
definite. We base our analysis upon the LaSalle invariance theorem [14, pp. 126-129].
A. Controller Design
Taking the time derivative of v (x) along the solutions of (2) results in
v̇ (x) = ∇v (x)T f (x) u
January 9, 2014

(9)
DRAFT

11

where ∇v (x) ∈ Rn . In order to get v̇ (x) ≤ 0, a direct way of designing u is proposed as
follows
u = −Q (x) f (x)T ∇v (x)

(10)

where Q : Rn → Rl×l and Q (x) ≥ ǫIl > 0, ǫ > 0, ∀x ∈ Rn . Then (9) becomes
v̇ (x) = −∇v (x)T f (x) Q (x) f (x)T ∇v (x) ≤ 0.

(11)

Substituting (10) into the continuous-time dynamical system (2) results in
ẋ = −f (x) Q (x) f (x)T ∇v (x)

(12)

with solutions which always satisfy the constraint c (x) = 0. The closed-loop system corresponding to the continuous-time dynamical system (2) and the controller (10) is depicted in
Fig.3.
Controller
x

u

Fig. 3.

Q x f x v x
T

System
x

f x u

Closed-loop control system

B. Convergence Analysis
Unlike a Lyapunov function, the objective function v (x) is not required to be positive
definite. As a consequence, the conclusions for Lyapunov functions are not applicable. Instead,
the invariance principle is applied to analyze the behavior of the solution of (12).
Theorem 4. Under Assumption 1, given x0 ∈ F , if the set K = {x ∈ Rn |v (x) ≤
v (x0 ) , c (x) = 0} is bounded, then the solution of (12) starting at x0 approaches x∗l ∈ S,
where S = {x ∈ K|∇v (x)T f (x) = 0}. If in addition V1 (x∗l ) = V2 (x∗l ) , then there must exist
P
∗
∗
∗
a λ∗ = [λ∗1 λ∗2 · · · λ∗m ]T ∈ Rm such that ∇v (x∗l ) = m
i=1 λi ∇ci (xl ) and c (xl ) = 0, namely

x∗l is a Karush–Kuhn–Tucker (KKT) point. Furthermore, if z T ∇xx L (x∗l , λ∗ ) z > 0, for all
P
z ∈ V1 (x∗l ) , z 6= 0, then x∗l is a strict local minimum, where L (x, λ) = v (x)− m
i=1 λi ci (x) .
Proof. The proof is composed of three propositions: Proposition 1 is to show that K is

compact and positively invariant with respect to (12); Proposition 2 is to show that the
January 9, 2014

DRAFT

12

solution of (12) starting at x0 approaches x∗l ∈ S; Proposition 3 is to show that x∗l ∈ S is
a KKT point, further a strict local minimum. The three propositions are proven in Appendix
E. 
Corollary 1. Suppose that f (x) is chosen as (7) for c : Rn → Rm and the set K =
{x ∈ Rn |v (x) ≤ v (x0 ) , c (x) = 0} is bounded for given x0 ∈ F . Then the solution of (12)
starting at x0 approaches x∗l ∈ S, where S = {x ∈ K|∇v (x)T f (x) = 0}, where x∗l is a
KKT point. In addition, if z T ∇xx L (x∗l , λ∗ ) z > 0, for all z ∈ V1 (x∗l ) , z 6= 0, then x∗l is a
P
strict local minimum, where L (x, λ) = v (x) − m
i=1 λi ci (x) .

Proof. Since V1 (x∗l ) = V2 (x∗l ) by Theorem 3, the remainder of the proof is the same as

that of Theorem 4. 
Corollary 2. Consider the following equality-constrained optimization problem
min v (x) , s.t. Ax = b.

x∈Rn

(13)

If (i) v (x) is convex and twice continuously differentiable, (ii) A ∈ Rp×n with rankA < n,
(iii) K ={x ∈ Rn |v (x) ≤ v (x0 ) , Ax = b} is bounded, then the solution of (12) with
f (x) = A⊥ starting at any x0 ∈ F approaches x∗ .
Proof. The solution of (12) starting at x0 approaches x∗l ∈ S. Since rankA < n, we have
V1 (x∗l ) = V2 (x∗l ) 6= ∅. Since the equality constrained optimization problem (13) is convex,
a KKT point x∗l is a global minimum x∗ of the problem (13). The remainder of proof is the
same as that of Theorem 4. 
Remark 6. If K is not a bounded set, then S defined in Theorem 4 may be empty.
Therefore, the boundedness of the set K is necessary. For example, v (x) = x1 + x2 , s.t.
c (x) = x1 − x2 = 0. The set K = { x ∈ R2 | x1 + x2 ≤ v (x0 ) , x1 − x2 = 0} is unbounded.
According to Theorem 1, we have f (x) = [1 1]T . In this case, ∇v (x)T f (x) = 2 6= 0 and
then the set S is empty.
C. A Modified Closed-Loop Dynamical System
Although the proposed approach ensures that the solutions satisfy the constraint, this
approach may fail if x0 ∈
/ F or if numerical algorithms are used to compute the solutions.
Moreover, if the impulse function δ is approximated, then the constraints will also be violated.
With these results, the following modified closed-loop dynamical system is proposed to amend
this situation.

January 9, 2014

DRAFT

13

Similar to [2], we introduce the term −ρ∇c (x) c (x) into (12), resulting in
ẋ = −ρ∇c (x) c (x) − f (x) Q (x) f (x)T ∇v (x) , x (0) = x0

(14)

where ρ > 0. Define vc (x) = c (x)T c (x) . Then
v̇c (x) = −ρc (x)T ∇c (x)T ∇c (x) c (x) ≤ 0,
where ∇c (x)T f (x) ≡ 0 is utilized. If the impulse function δ is approximated, then ∇c (x)T f (x) ≈
0 and can be ignored in practice. Therefore, the solutions of (14) will tend to the feasible
set F if ∇c (x) is of full column rank. Once c (x) = 0, the modified dynamical system
(14) degenerates to (12). The self-correcting feature enables the step size to be automatically
controlled in the numerical integration process or to tolerate uncertainties when the differential
equation is realized by using analog hardware.
Remark 7. The matrix Q (x) plays a role in coordinating the convergence rate of all
states by minimizing the condition number of the matrix functions like f (x) Q (x) f (x)T .
Moreover, it also plays a role in avoiding instability in the numerical solution of differential
equations by normalizing the Lipschitz condition of functions like f (x) Q (x) f (x)T ∇v (x) .
Concrete examples are given in the following section.
V. I LLUSTRATIVE E XAMPLES
A. Estimate of Attraction Domain
For a given Lyapunov function, the crucial step in any procedure for estimating the
attraction domain is determining the optimal estimate. Consider the system of differential
equations:
ẋ = Ax + g (x)

(15)

where x ∈ Rn is the state vector, A ∈ Rn×n is a Hurwitz matrix, and g : Rn → Rn is a
vector function. Let v (x) = xT P x be a given quadratic Lyapunov function for the origin of
(15), i.e., P ∈ Rn×n is a positive-definite matrix such that AT P + P A < 0. Then the largest
ellipsoidal estimate of the attraction domain of the origin can be computed via the following
equality-constrained optimization problem [15]:
min xT P x s.t. xT P [Ax + g (x)] = 0.

x∈Rn \{0}

Since {x ∈ Rn |xT P x ≤ xT0 P x0 } is bounded, the subset
K = {x ∈ Rn |xT P x ≤ xT0 P x0 , xT P [Ax + g (x)] = 0}
January 9, 2014

DRAFT

14

is bounded no matter what g is.
For simplicity, consider (15) with x = [x1 x2 ]T ∈ R2 , A = −I2 , P = I2 and g (x) =

(σ (x) + 1) [x1 x2 ]T , where σ (x) = (x1 + x2 + 2) (x2 + 1) − 0.1 (x1 + 1)2 . Then the optimization problem is formulated as


2
2
2
2
min
x
+
x
s.t.
x
+
x
σ (x) = 0.
1
2
1
2
2

x∈R \{0}

Since x 6= 0, the problem is further formulated as

min2 v (x) = x21 + x22 s.t. σ (x) = 0.

x∈R

Then
∇v (x) = [2x1 2x2 ]T


d2 − 0.1d21 − 0.2d1 d3

∇c (x) = 
d2 − 0.1d21 + d3

d1 = x1 + 1, d2 = x2 + 1, d3 = x1 + x2 + 2.

In this example, we adopt the modified dynamics (14), where f is chosen as (7) with δ (x) =

e−γ|x| , and the parameters are chosen as γ = 10, ρ = Q = 20 ∇cc − f f T ∇v . We solve
the differential equation (14) by using the MATLAB function “ode45” with “variable-step1”.

Compared with the MATLAB optimal constrained nonlinear multivariate function “fmincon”,
we derive the comparisons in Table 1.
The point xs = [−1 −1]T is a singular point, at which ∇c (xs ) = [0 0]T . As shown
in Table 1, under initial points [−3 1]T ∈ F and [2 −4]T ∈ F , the MATLAB function
fails to find the minimum and stops at the singular point, whereas the proposed approach
still finds the minimum. Under initial point [1 −4]T ∈
/ F , the proposed approach can still
find the minimum, similar to the MATLAB function. Under a different initial value, the
evolutions of (14) are shown in Fig.4. As shown, once close to the singular point [−1 −1]T , the
solutions of (14) change direction and then move to the minimum x∗l = [0.2061 − 0.8545]T .
Compared with the discrete optimal methods offered by MATLAB, these results show that the
proposed approach avoids convergence to a singular point. Moreover, the proposed approach
is comparable with currently available conventional discrete optimal methods and facilitates
even faster convergence. The latter conclusion is consistent with that proposed in [1],[3].
1

In this section, all computation is performed by MATLAB 6.5 on a personal computer (Asus x8ai) with Intel core Duo

2 Processor at 2.2GHz.
January 9, 2014

DRAFT

15

TABLE 1. COMPUTED RESULT FOR EXAMPLE 1
Method

Initial Point

Solution

Optimal Value

cpu time (sec.)

Matlab fmincon

[-3 1]T

[-1 -1]T

2.0000

Not Available

New method

[-3 1]T

[0.2062 -0.8546]T

0.7729

0.125

Matlab fmincon

[2 -4]T

[-1 -1]T

2.0000

Not Available

New method

[2 -4]T

[0.2062 -0.8545]T

0.7726

0.0940

Matlab fmincon

[1 -4]T

[0.2143 -0.8533]T

0.7740

0.2030

New method

[1 -4]T

[0.2056 -0.8550]T

0.7733

0.1100

.

1
Evolution trajectory of
initial point (−3,1)
0.5

0
(0,0)
−0.5

−1

(0.2061 −0.8545)

x2

Singular point (−1,−1)
−1.5

−2
Evolution trajectory of
initial point (2,−4)

−2.5

−3
Evolution trajectory of
initial point (1,−4)

−3.5

−4
−3

Fig. 4.

−2

−1

0
x1

1

2

3

Optimization for estimate of attraction domain. Solution Evolution (solid line), Constraint (dot line), Objective

(dash-dot line).

January 9, 2014

DRAFT

16

B. Estimate of Essential Matrix
For simplicity, assume that images are taken by two identical pin-hole cameras with focal
length equal to one. The two cameras are specified by the camera centers C1 , C2 ∈ R3 and
attached orthogonal camera frames {e1 , e2 , e3 } and {e′1 , e′2 , e′3 }, respectively. Denote T =
C2 − C1 ∈ R3 to be the translation from the first camera to the second and R ∈ R3×3 to be
the rotation matrix from the basis vectors {e1 , e2 , e3 } to {e′1 , e′2 , e′3 }, expressed with respect
to the basis {e1 , e2 , e3 } . Then, it is well known in the computer vision literature [16] that
two corresponding image points are represented as follows:
1
Mk ,
Mk (3)
1
= ′
Mk′ , k = 1, 2, · · · , N
Mk (3)

m1,k =
m2,k

(16)

where Mk , Mk′ represent the positions of the kth point expressed in the two camera frames
{e1 , e2 , e3 } to {e′1 , e′2 , e′3 } , respectively; Mk (3) , Mk′ (3) represent the third element of vectors
Mk , Mk′ , respectively. They have the relationship Mk = RMk′ + T, k = 1, 2, · · · , N. These
corresponding image points satisfy the socalled epipolar constraint [16, p. 257]:
mT1,k Em2,k = 0, k = 1, 2, · · · , N

(17)

where E = [T ]× R is known as the essential matrix.
Mk

m2,k

m1,k

e1

C1

e1a

e2

e2a
C2

e3

e3a

R, T

Fig. 5.

Epipolar geometry

By using the direct product ⊗ and the vec (·) operation, the equations in (17) are equivalent
to
Aϕ = 0N ×1

January 9, 2014

(18)

DRAFT

17

where


mT2,1



A=


⊗
..
.

mT1,1

mT2,N ⊗ mT1,N

ϕ = vec [T ]× R .





 ∈ RN ×9 ,


(19)

In practice, these image points m1,k and m2,k are subject to noise, k = 1, 2, · · · , N. Therefore,
T and R are often solved by the following optimization problem
1
min12 v (x) = ϕ (x)T AT Aϕ (x)
x∈R
2

1
s.t.
kT k2 − 1 = 0
2

1 T
R R − I3 = 03×3
2

(20)

where x = [T T vecT (R)]T ∈ R12 . This is an equality-constrained optimization considered
here. In the following, the proposed approach is applied to the optimization problem (20).

By Theorem 2, the projection matrix for the constraint 12 kT k2 − 1 = 0 is
f = I3 −

TTT

.
δ kT k2 + kT k2


Since kT k2 = 1 has to be satisfied exactly or approximately, then δ kT k2 = 0. So, the
projection matrix for the constraint is


f = I3 − T T kT k2 .

Then the constraint is transformed into



Ṫ = I3 − T T T kT k2 u1,

where u1 ∈ R3 . By (3), the constraint

1
2


RT R − I3 = 03×3 is transformed into

Ṙ = [u2 ]× R,

where u2 ∈ R3 . Furthermore, the equation above is rewritten as
 

vec Ṙ = RT ⊗ I3 Hu2 .

January 9, 2014

DRAFT

18

Then the continuous-time dynamical system, whose solutions always satisfy the equality


constraints 12 kT k2 − 1 = 0 and 21 RT R − I3 = 03×3 , is expressed as (2) with



I3 − T T T kT k2
03×3
12×6
f (x) = 
,
 ∈R
T
09×3
R ⊗ I3 H


u1
 ∈ R6 .
u=
(21)
u2

If the initial value kT (0)k2 = 1 and R (0)T R (0) = I3 , then all solutions of (2) satisfy the

equality constraints. Since ∇v (x) = [ RT ⊗ I3 H I3 ⊗ [T ]× ]T AT Aϕ, the time derivative
of v (x) along the solutions of (2) is

v̇ (x) = −ϕT AT AΘ (x)T Q (x) Θ (x) AT Aϕ ≤ 0,
where




T
T 
I3 − T T T kT k2 H T RT ⊗ I3
 ∈ R6×9 .
Θ (x) = 
T
T
T
T
H R ⊗ I3
I3 ⊗ [T ]×

The simplest way of choosing Q (x) is Q (x) ≡ I6 . In this case, the eigenvalues of the
matrix AΘT (x) Θ (x) AT are often ill-conditioned, namely


λmin AΘT (x) Θ (x) AT ≪ λmax AΘT (x) Θ (x) AT .

Convergence rates of the components of Aϕ (x) depend on the eigenvalues of AΘT (x) Q (x) Θ (x) AT .
As a consequence, some components of Aϕ converge fast, while the other may converge
slowly. This leads to poor asymptotic performance of the closed-loop system. It is expected
that each component of Aϕ can converge at the same speed as far as possible. Suppose that
there exists a Q̄ (x) such that
AΘT (x) Q̄ (x) Θ (x) AT = I9 .
Then
v̇ (x) ≤ −ϕT AT Aϕ ≤ 0.
By Theorem 4, x will approach the set { x ∈ Rn | Aϕ (x) = 0} , each element of which is a
global minimum since v (x) = 0 in the set. Moreover, each component of Aϕ converges at
a similar speed. However, it is difficult to obtain such a Q̄ (x), since the number of degrees
of freedom of Q̄ (x) ∈ R6×6 is less than the number of elements of I9 . A modified way is
to make AΘT (x) Q (x) Θ (x) AT ≈ I9 . A natural choice is proposed as follows


†
T
T
Q (x) = µ Θ (x) A AΘ (x)
+ ǫI6
January 9, 2014

(22)
DRAFT

19

where µ > 0, Θ (x) AT AΘT (x)

†

denotes the Moore Penrose inverse of Θ (x) AT AΘT (x).

The matrix ǫI6 is to make Q (x) positive definite, where ǫ is a small positive real. From the
†
procedure above, Θ (x) AT AΘT (x) needs to be computed every time. This however will
cost much time. A time-saving way is to update Q (x) at a reasonable interval. Then (12)

becomes


†
T
T
ẋ = −µf (x)
Θ (x) A AΘ (x)
+ ǫI6 Θ (x) AT Aϕ (x)

(23)

where f (x) is defined in (21). The differential equation can be solved by Runge-Kutta
methods, etc. The solutions of (23) satisfy the constraints, where x = [T T vec(R)T ]T .
Moreover, the dynamic system will reach some final resting state eventually.
Suppose that there exist 6 points in the field of view, whose positions are expressed in the
first camera frame as follows: M1 = [−1 1 1]T , M2 = [2 0 1]T , M3 = [1 −1 1]T , M4 = [−1
−1 1]T , M5 = [1 1 1]T , M6 = [−1 3 1]T . Compared with the first camera frame, the second
camera frame has translated and rotated with




1
0.9900
−0.0894
0.1088








T̄ =  1  , R̄ =  0.0993
0.9910 −0.0894  .




−1
−0.0998 0.0993
0.9900

The image points are generated by (16). Using the generated image points, we obtain A by

(19). Setting the initial value as follows T (0) = [0 0 1]T , R (0) = I3 , µ = 20, ǫ = 0.01. We
solve the differential equation (14) by using MATLAB function “ode45” with “variable-step”.
Compared with MATLAB optimal constrained nonlinear multivariate function “fmincon”, we
have the following comparisons:
TABLE 2. COMPUTED RESULT FOR EXAMPLE 2
Method

R∗T R̄ − I3

cpu time (sec.)

MATLAB fmincon

1.2469e-004

0.2500

New Approach

1.8784e-005

0.1400

.

As shown in Table 2, the proposed approach requires less time to achieve a higher accuracy.
Given that v (x∗ ) = 0, the solution is a global minimum. The evolution of each element
of x is shown in Fig.5. The state eventually reaches a rest state at a similar speed. With
different initial values, several other simulations are also implemented. Based on the results,
the proposed algorithm has met the expectations.

January 9, 2014

DRAFT

20

1

0.8

0.6

0.4

0.2

0

−0.2

−0.4

−0.6

0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

0.45

0.5

time (sec)

Fig. 6.

Evolvement of the state

VI. C ONCLUSIONS
An approach to continuous-time, equality-constrained optimization based on a new projection matrix is proposed for the determination of local minima. With the transformation of the
equality constraint into a continuous-time dynamical system, the class of equality-constrained
optimization is formulated as a control problem. The resultant approach is more general than
the existing control theoretic approaches. Thus, the proposed approach serves as a potential
bridge between the optimization and control theories. Compared with other standard discretetime methods, the proposed approach avoids convergence to a singular point and facilitates
faster convergence through numerical integration on a digital computer.
A PPENDIX
A. Kronecker Product and Vec
The symbol vec(X) is the column vector obtained by stacking the second column of X
under the first, and then the third, and so on. With X = [xij ] ∈ Rn×m , the Kronecker product
X ⊗ Y is the matrix





X⊗Y =


x11 Y
..
.
xn1 Y

···
..
.



x1m Y


..
.
.

· · · xnm Y


In fact, we have the following relationships vec(XY Z) = Z T ⊗ X vec(Y ) [17, p. 318].
January 9, 2014

DRAFT

21

B. Skew-Symmetric Matrix
The cross product of two vectors x ∈ R3 and y ∈ R3 is denoted by x × y = [x]× y, where
the symbol [·]× : R3 → R3×3 is defined as [13, p. 194]:


0
−x
x
3
2




3×3
[x]× ,  x3
0 −x1  ∈ R .


−x2 x1
0

By the definition of [x]× , we have x × x = [x]× x = 03×1 , ∀x ∈ R3 and

vec [x]× = Hx,

T
0 0 1 0 −1 0 
 0 0 0


H =  0 0 −1 0 0 0 1 0 0  .


0 1 0 −1 0 0 0 0 0

C. Proof of Lemma 1


Since δ kLk2 + kLk2 = 1 if L = 0 and δ kLk2 + kLk2 = kLk2 if L 6= 0, we have

δ kLk2 + kLk2 6= 0, ∀L ∈ Rn . According to this, we have the following relationship


δ kLk2 + kLk2



= LT − LT kLk2 δ kLk2 + kLk2

LT In − LLT



≡ 0, ∀L ∈ Rn .

This implies that LT z = 0, ∀z ∈ W2 , namely W2 ⊆ W1 . On the other hand, any z ∈ W1 is
rewritten as
z = In − LLT





δ kLk2 + kLk2 z

where LT z = 0 is utilized. Hence W1 ⊆ W2 . Consequently, W1 = W2 .
D. Proof of Theorem 3
Denote
V1j = {z ∈ Rn |∇cTi z = 0, i = 1, · · · , j, j ≤ m}
V2j = {z ∈ Rn |z = fj uj , uj ∈ Rn , j ≤ m}.
First, by Theorem 2, it is easy to see that the conclusions are satisfied with j = 1. Assume
V1k−1 = V2k−1 and then prove that V1k = V2k holds. If so, then we can conclude this proof.

January 9, 2014

DRAFT

22

By V1k−1 (x) = V2k−1 (x) , we have
V1k = {z ∈ Rn |∇cTk z = 0, z ∈ V1k−1 }
= {z ∈ Rn |∇cTk z = 0, z = fk−1 uk−1, uk−1 ∈ Rn }
= {z ∈ Rn |∇cTk fk−1 uk−1 = 0, z = fk−1 uk−1 , uk−1 ∈ Rn }.
By Lemma 1, we have
∇cTk fk−1 uk−1 = 0 ⇔

uk−1

namely,

T
fk−1
∇ck ∇cTk fk−1



= In −
2
T
T
∇ck
+ fk−1
δ fk−1 ∇ck

2



 uk ,

V1k = V2k = {z ∈ Rn |z = fk uk , uk ∈ Rn }


T ∇c ∇cT f
fk−1
k−1
k
k

where fk = fk−1 In −  T
.
2
2
T ∇c
δ kfk−1 ∇ck k +kfk−1
kk
E. Proof of Propositions in Theorem 3
(i) Proof of Proposition 1. In the space Rn , the set K is compact iff it is bounded and closed
by Theorem 8.2 in [18, p.41]. Hence, the remainder of work is to prove that K is closed.
Suppose, to the contrary, K is not closed. Then there exists a sequence x (tn ) ∈ K → p ∈
/K
with tn → ∞. Whereas, v (p) = lim v (x (tn )) ≤ v (x0 ) and c (p) = lim c (x (tn )) = 0
tn →∞

tn →∞

which imply p ∈ K. The contradiction implies that K is closed. Hence, the set K is compact.
By (11), v (x) ≤ v (x0 ) with respect to (12), t ≥ 0. By Assumption 1, all solutions of (12)
satisfy c (x) = 0. Therefore, K is positively invariant with respect to (12).
(ii) Proof of Proposition 2. Since K is compact and positively invariant with respect to
(12), by Theorem 4.4 (invariance principle) in [14, p. 128], the solution of (12) starting at
x0 approaches v̇ (x) = 0, namely ∇v (x)T f (x) = 0. In addition, since (12) becomes ẋ = 0
in S, the solution approaches a constant vector x∗l ∈ S.
(iii) Proof of Proposition 3. Since V1 (x∗l ) = V2 (x∗l ) and x∗l ∈ S satisfy the following two
equalities
∇v (x∗l )T f (x∗l ) = 0, c (x∗l ) = 0,
there exists a u such that z = f (x∗l ) u for any z ∈ V1 (x∗l ) . As a consequence, for any
z ∈ V1 (x∗l ) , ∇v (x∗l )T z = ∇v (x∗l )T f (x∗l ) u = 0. There must exist λ∗i ∈ R, i = 1, · · · , m
P
∗
∗
∗
∗ T
such that ∇v (x∗l ) = m
i=1 λi ∇ci (xl ). Otherwise ∃z̄ ∈ V1 (xl ), ∇v (xl ) z̄ 6= 0. Therefore,
January 9, 2014

DRAFT

23

x∗l ∈ S is a KKT point [12, p.328]. Furthermore, by Theorem 12.6 in [12, p.345], x∗l is a
strict local minimum if z T ∇xx L (x∗l , λ∗ ) z > 0, for all z ∈ V1 (x∗l ) , z 6= 0.

R EFERENCES
[1] K. Tanabe, A geometric method in nonlinear programming, Journal of Optimization Theory and Applications. 30(1980)
181–210.
[2] H. Yamashita, A differential equation approach to nonlinear programming, Mathematical Programming, 18 (1980),
155–168.
[3] A.A. Brown, M.C. Bartholomew-Biggs, ODE versus SQP methods for constrained optimization, Journal of Optimization Theory and Applications, 62 (1989), 371–386.
[4] S. Zhang, A.G. Constantinides, Lagrange programming neural networks, IEEE Transactions on Circuits and Systems-II:
Analog and Digital Signal Processing, 39 (1992), 441–452.
[5] Z.-G. Hou, A hierarchical optimization neural network for large-scale dynamic systems, Automatica, 37 (2001), 1931–
1940.
[6] L.-Z. Liao, H. Qi, L. Qi, Neurodynamical optimization, Journal of Global Optimization, 28 (2004), 175–195.
[7] M.P. Barbarosou, N.G. Maratos, A nonfeasible gradient projection recurrent neural network for equality-constrained
optimization problems, IEEE Transactions on Neural Networks, 19 (2008), 1665–1677.
[8] P.-A. Absi, Computation with continuous-time dynamical systems, in the Grand Challenge in Non-Classical Computation International Workshop, York, United Kingdom, 2005, Apr. 18–19.
[9] J.J. Hopfield, D.W. Tank, Neural computation of decisions in optimization problems, Biological Cybernetics, 52 (1985),
141–152.
[10] D.G. Luenberger, Y. Ye, Linear and Nonlinear Programming, third ed., Springer, Boston, 2008.
[11] E.L. Lawler, D.E. Wood, Branch-and-bound methods: a survey, Operations Research, 14 (1966), 699–719.
[12] J. Nocedal, S.J. Wright, Numerical Optimization, Springer-Verlag, New York, 1999.
[13] A. Isidori, L. Marconi, A. Serrani, Robust Autonomous Guidance: An Internal Model-Based Approach, SpringerVerlag, London, 2003.
[14] H.K. Khalil, Nonlinear Systems, third ed., Prentice-Hall, Upper Saddle River, New York, 2002.
[15] G. Chesi, A. Garulli, A. Tesi, A. Vicino, Solving quadratic distance problems: an LMI-Based approach, IEEE
Transaction on Automatic Control, 48 (2003), 200–212.
[16] R. Hartley, A. Zisserman, Multiple View Geometry in Computer Vision, second ed., Cambridge University Press,
Cambridge, 2003.
[17] U. Helmke, J.B. Moore, Optimization and Dynamical Systems. Springer-Verlag, 1994.
[18] F. Morgan, Real Analysis and Applications: Including Fourier Series and the Calculus of Variations. American
Mathematical Society, 2005.

January 9, 2014

DRAFT

x

f x f x

T

x

v x

c x

0

v x

v x

r

