Quadratic Programming forMonotone Control Theoretic Splines
Masaaki Nagahara

Clyde F. Martin

Yutaka Yamamoto

Graduate School of Informatics
Kyoto University
Kyoto, Japan
Email: nagahara@ieee.org

Department of Mathematics & Statistics
Texas Tech University
Texas, USA
Email: clyde.f.martin@ttu.edu

Graduate School of Informatics
Kyoto University
Kyoto, Japan
Email: yy@i.kyoto-u.ac.jp

Abstract—In this article, we solve the problem of monotone
control theoretic splines. Monotone control theoretic splines havebeen solved only when the target system is the second-orderintegrator 1/s2 , but not for other cases. To solve this, we first
formulate the problem as semi-infinite quadratic programming,and then we adopt discretization technique. The problem isreduced to a finite-dimensional quadratic programming, which
can be easily solved by numerical computation. We also illustrateexamples to show the effectiveness of our method.

I. I NTRODUCTION
Control theoretic splines are interpolating smoothing splines
with a constraint written in linear differential equations [4],that is, the interpolating curve between data points is obtainedby an output of a given linear system. They have been proved
to be useful in trajectory planning for aircrafts and robots, andstatistics to smooth noisy data [3].
On the other hand, monotone control theoretic splines arealso important in deriving a model or an estimation of aparameter such as the growing rate of an individual [2].
This needs to add a non-negative derivative constraint. The
difficulty of this problem is that this constraint is infinite
dimensional. In [2], this has been solved when the target linear
system is a second-order integrator 1/s2 . However, for generalcases, there has been no solution.
In this article, we propose a new method for monotone
splines with general linear systems. To solve this, we first
formulate the problem as a semi-infinite quadratic programming. Then we adopt discretization technique [5] to reduce theinfinite-dimensional constraints to a finite dimensional ones.
This discretization technique is also know as fast discretization in sampled-data control theory [1]. We will show that
our approximation satisfies the original infinite-dimensional
constraints provided that the number of discretization is sufficiently large. The finite-dimensional quadratic programming
can be easily solved by numerical computation softwares,for example, MATLAB Optimization Toolbox [6]. We also
illustrate an example to show the effectiveness of our method.
Notations
In this article, we use the following notations. R, Rn , andare respectively the set of real numbers, n-dimensionalR
real vectors and n × n matrices. We denote L2 [0, T ] byn×n

the Lebesgue space consisting of all square integrable realfunctions on [0, T ] ⊂ R, endowed with the inner product(x, y) := Tx, y ∈ L2 [0, T ].

x(t)y(t)dt,0

For a matrix (a vector) M , M  is the transpose of M . Fora vector v = [v1 , . . . , vn ] ∈ Rn , v  0 means vi ≤ 0 for
i = 1, 2, . . . , n, and for two vectors v and w, v  w iffv − w  0. For a vector v ∈ Rn ,v1 :=n|vi |.

i=1

II. M ONOTONE C ONTROL T HEORETIC S PLINES
We consider the following problem of monotone control
theoretic splines:Problem 1 (Monotone control theoretic splines): Given alinear time-invariant systemẋ = Ax + Bu,

y = Cx,x(0) = 0,

where A ∈ Rn×n , B ∈ Rn×1 and C ∈ R1×n , and also given
a data set {(t1 , α1 ), (t2 , α1 ), . . . (tN , αN )}, with 0 < t1 <t2 < · · · < tN =: T , and α1 , . . . , αN ∈ R, find the controlu ∈ L2 [0, T ] that minimizes the following cost function:J(u) := λ T02u(t) dt +Nwi (y(ti ) − αi )2 ,(1)i=1where wi ’s are given positive numbers (weights), with themonotonicity constraintẏ(t) ≥ 0,∀t ∈ [0, T ].

(2)This problem has been solved in [2] only when 
0 10
A=, B=, C= 1 0 .
0 0
1However, this problem has not been solved for other cases.

III. F ORMULATION BY SEMI - INFINITE QUADRATICPROGRAMMING

We here formulate our problem by semi-infinite quadratic
programming. We first assume that we search the optimalcontrol u in the following subspace of L2 [0, T ] as taken in[4]:N
θi gti , θi ∈ R ,VN := u ∈ L2 [0, T ] : u =The constraint (2) now becomes Φ(t) θ ≥ 0 for all t ∈ [0, T ].
In summary, Problem 1 is formulated by the following.
Problem 2: Find θ ∈ RN that minimize
1f (θ) = θ P θ + q  θ,
2
P := 2 (λI + GW ) G, q = −2W Gα,such thatg(θ, t) = −Φ(t) θ ≤ 0,i=1

where gt is defined byt > τ,t ≤ τ.

IV. D ISCRETIZATION APPROACH TO MONOTONE
THEORETIC CONTROL SPLINES

We also denote the inner product in L2 [0, T ] by (·, ·). By
assuming that the control u ∈ VN , the cost function (1)becomesNJ
θi gti = θ (λI + GW ) Gθ − 2αGW θ + αW α,i=1

where α := [α1 , . . . , αN ] , θ := [θ1 , . . . , θN ] , and⎤⎡(gt1 , gt1 ) (gt2 , gt1 ) . . . (gtN , gt1 )
⎥
⎢
..
..
..
..
G := ⎣⎦,.
.
.
.
(gt1 , gtN ) (gt2 , gtN ) . . .
W := diag {w1 , . . . , wN } .

=N

i=1eA(t−τ ) B

N

i=1θi

A(t−τ )

 t

CAe0g(θ, Ti ) ≤ −,N

θi gti (t)i=1

Bgti (τ )dτ + CBgti (t) .

By using this, ẏ(t) can be calculated as

NT

θi
ġt (τ )gti (τ )dτ + CBgti (t)ẏ(t) =0N=θi (ġt , gti ) + CBgti (t)i=1=: Φ(t) θ,where⎡⎤⎡i = 0, 1, 2, . . . , M,−λ1  θ  λ1,Let ġt (τ ) be the derivative of gt (τ ), that is,
CAeA(t−τ ) B, if t > τ∂gt (τ )=ġt (τ ) :=∂t0, if t ≤ τ,i=1

where Ti ’s are discretization points satisfying
Then we evaluate the function g(θ, t) in (3) at the discretization points t = T0 , T1 , . . . , TM . To guarantee the constraint(3), we choose real numbers  > 0 and λ > 0 and adopt thefollowing constraints:

θi gti (τ )dτ + CB

0[0, T ] = [T0 , T1 ] ∪ (T1 , T2 ] ∪ · · · ∪ (TM −1 , TM ],0 = T0 < T1 < · · · < TM −1 < TM = T.

ẏ(t) = C ẋ(t)= CAx(t) + CBu(t) t= CAeA(t−τ ) Bu(τ )dτ + CBu(t)= CA

The difficulty of Problem 2 is that the inequality constraint(3) is infinite. We here introduce an approximation technique
for such an infinite-dimensional constraint. To reduce this to a
finite dimensional one, we adopt the technique of discretization[5], which is also known as fast discretization [1] in sampleddata control theory.
We divide the interval [0, T ] into M subintervals:(gtN , gtN )Then the derivative ẏ(t) is calculated as0(3)This is an semi-infinite quadratic programming.

CeA(t−τ ) B,gt (τ ) :=0, t

for all t ∈ [0, T ].

⎤gt1 (t)(ġt , gt1 )⎢ gt2 (t) ⎥⎢ (ġt , gt2 ) ⎥⎢⎢⎥⎥Φ(t) := ⎢⎥ + CB ⎢ .. ⎥ .
..
⎣ . ⎦⎣⎦.
(ġt , gtN )gtN (t).

(4)

where 1 = [1, 1, . . . , 1] ∈ RN and  denotes a component
wise inequality, that is, the second inequality means −λ ≤θi ≤ λ for all i = 1, 2, . . . , N . Define feasible sets F andF,M by
F := θ ∈ RN : θ satisfies (3). ,
F,M := θ ∈ RN : θ satisfies (4). .
Let
Imax (M ) :=max

i=1,2,...,M|Ti − Ti−1 |,μ := max Φ(t)1 .
t∈[0,T ]Then we have the following proposition:
Proposition 1: Assume 0 < μ < ∞. Then F,M ⊆ F if.
(5)
Imax (M ) ≤λμProof. Let θ ∈ F,M . By (4), θ satisfiesg(θ, Ti ) ≤ −,i = 0, 1, . . . , M.

Then, for any t ∈ [0, T ], there exists i ∈ {1, . . . , M } such thatt ∈ [Ti−1 , Ti ]. Since ∂g(θ, t) =Φ(t) θ ≤ λμ.
maxmax∂tt∈[0,T ], θ∈Ft∈[0,T ], θ∈Fwe have

VI. C OMPUTING I NNER P RODUCTg(θ, t) − g(θ, Ti ) ≤ λμ(Ti − t).

By this inequality, we haveg(θ, t) ≤ λμ(Ti − t) + g(θ, Ti )≤ λμ(Ti − t) − ≤ λμImax (M ) − ≤ 0.
and hence, θ ∈ F. Therefore F,M ⊆ F.

By this proposition, we can guarantee the constraint (3) bysearching the parameter θ in the finite-dimensional feasible
set F,M provided that the number M is sufficiently large tosatisfy (5).
The finite-dimensional constraints (3) can be represented asa matrix inequality,Hθ  v,
where⎡ ⎤−⎢−⎥⎢ ⎥⎢ .. ⎥⎢ ⎥v := ⎢ . ⎥ .
⎢−⎥⎢ ⎥⎣λ1⎦λ1⎡⎤Φ(T0 )
⎢ Φ(T1 ) ⎥⎢⎥⎢⎥..
⎢⎥.
H := ⎢⎥,⎢Φ(TM ) ⎥⎢⎥⎣⎦I−IIn summary, the problem of the monotone control theoretic spline for arbitrarily given linear time-invariant system{A, B, C} (note that no assumption is needed on A, B andC) can be solved by the following quadratic programming:
find θ ∈ RN which minimizesf (θ) =1 θ P θ + q  θ,
2subject toHθ  v.
This is a standard quadratic programming and can be easily
solved by MATLAB (using the command quadprog inMATLAB Optimization Toolbox [6]).
V. E QUALITY CONSTRAINTS
We can also include equality constraints into our optimization. For example, if we want to have y(T ) = 1 and ẏ(T ) = 0,then our constraints are represented by⎡⎤⎤⎡(gT , gt1 )(ġT , gt1 )
⎢
⎥
⎥
⎢
..
..
⎣⎦ θ = 1, ⎣⎦ θ = 0.
.
.
(gT , gTN )(ġT , gTN )These constraints are finite dimensional and can be included
easily in our quadratic programming. In general, equality
constraints on y(t) and ẏ(t), t ∈ [0, T ] can be represented
linear constraints as above.

To compute the Grammian G and the matrix Φ(t), we haveto compute the following inner products for fixed s, t ∈ [0, T ]:
 T
 T
(gs , gt ) =gs (τ )gt (τ )dτ, (ġs , gt ) =ġs (τ )gt (τ )dτ.
00

These values can be easily computed by matrix exponentials[1, Lemma 10.5.1]:
 T
gs (τ )gt (τ )dτ(gs , gt ) =0
 min(s,t)  
CeA(s−τ ) B
CeA(t−τ ) B dτ=0= B  eA s min(s,t)
0e−A τ C  Ce−Aτ dτeAt B 
F12 vt ,=: vs F22where vτ := eAτ B (τ = t, s), and the matrices F22 and F12are defined by  
F11 F12AC C:= exp
h , h := min(s, t).
0 F220−A
The inner product (ġs , gt ) is also obtained by (ġs , gt ) = vs H22H12 vt ,  H11 H12AA C  C:= exph .
0H220−AVII. E XAMPLESWe here show an example of monotone control theoretic
spline. Assume the system is given by⎡⎤⎡ ⎤0 1 00
A = ⎣0 0 1⎦ , B = ⎣0⎦ , C = 1 0 0 .
0 0 01This system has the transfer function 1/s3 . Let the original
curve be yorig (t) = 1 − e−t , and the data is obtained bysampling yorig (t) to which Gaussian noise with N (0, 0.01) isindependently added. The estimated parameter is given byθ = [93.9, 185.1, 17.1, − 24.4, − 7.2, 12.1,− 10.5, 9.7, − 4.7, 0.4, 0.3]

(6)

The estimation result is shown in Fig. 1. This figure shows that
our estimation works well. The derivative of the estimationy(t) is given in Fig. 2. We can see that the derivative of y(t)is always non-negative for all t ∈ [0, T ].
We then consider another system given by 

1 01
A=, B=, C= 0 1 .
1 00Note that the transfer function of this system is 1/s(s−1). For
this system, we obtain the optimal θ by the proposed method.
The parameter θ is obtained asθ = [15.05, 9.20, 0.52, − 1.70, 1.07, − 0.52,− 0.05, 0.20, 0.24, − 0.40, 0.12, ] .

0.25
10.9
0.2
0.8
0.7
0.15

0.6
0.5

0.1

0.4
0.3

0.05

0.2
0.1
00123

45006Fig. 1. Data (circles), original curve yorig (t) = 1−e−t (dash), and estimationy(t) (solid) for the system 1/s3Fig. 4.

12

3

4

5

6

Absolute estimation errors for 1/s3 (dash) and 1/s(s − 1) (solid)

Note that the values of the elements of this θ is much smaller
than that in (6). The estimation result is shown in Fig. 3. We
can see that the estimation for the system 1/s(s − 1) is better
than that for 1/s3 . This is because the system 1/s(s−1) has afaster mode s = 1. Fig. 4 shows the absolute estimation errors,that is, |yorig (t)−y(t)|. By this figure, the error for 1/s(s−1) ismuch reduced compared with that for 1/s3 . This example also
shows that for given data, the choice of the system {A, B, C}
considerably affects the accuracy of estimation. How to choose
the best {A, B, C} for given data is a future work of this study.

derivative of y
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2

VIII. C ONCLUSION

0.1

In this article, we have proposed a new method for solving
the problem of monotone control theoretic splines for generallinear systems. By using fast discretization technique, theproblem is described as quadratic programming, and hence
the optimal parameters can be easily obtained. Splines withanother constraint such as concavity, i.e., ÿ(t) ≥ 0 can be also
solved by the same method.

00Fig. 2.

123456

Derivative of estimation y(t) for the system 1/s3

R EFERENCES
1[1] T. Chen and B. Francis, Optimal Sampled-Data Control Systems,Springer, 1995.
[2] M. Egerstedt and C.F. Martin, Monotone smoothing splines, Proceedingsof the Mathematical Theory of Networks and Systems, 2000.
[3] M. Egerstedt and C.F. Martin, Control Theoretic Splines, Princeton
University Press, 2010.
[4] S. Sun, M. B. Egerstedt, and C. F. Martin, Control theoretic smoothingsplines, IEEE Trans. on Automatics Control, Vol. 45, No. 12, pp. 2271–2279, 2000.
[5] K. L. Teo and X. Q. Yang, Computational discretization algorithms
for functional inequality constrained optimization, Annals of Operations
Research, Vol. 98, pp. 215–234, 2000.
[6] http://www.mathworks.com/products/optimization/.

0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
0

1

2

3

4

56Fig. 3. Data (circles), original curve yorig (t) = 1−e−t (dash), and estimationy(t) (solid) for the system 1/s(s − 1)