Abstract
A computational method based on Bézier control points is presented to solve optimal control problems governed by time varying linear dynamical systems subject to terminal state equality constraints and state inequality constraints. The method approximates each of the system state variables and each of the control variables by a Bézier curve of unknown control points. The new approximated problems converted to a quadratic programming problem which can be solved more easily than the original problem. Some examples are given to verify the efficiency and reliability of the proposed method. Mathematical subject classification: 49N10.
Bézier control points; constrained optimal control; linear time varying dynamical systems
Bézier control points method to solve constrained quadratic optimal control of time varying linear systems
F. Ghomanjani; M.H. Farahi; M. Gachpazan
Department of Applied Mathematics, Faculty of Mathematical Sciences, Ferdowsi University of Mashhad, Mashhad, Iran. Email: gachpazan@ferdowsi.um.ac.ir
ABSTRACT
A computational method based on Bézier control points is presented to solve optimal control problems governed by time varying linear dynamical systems subject to terminal state equality constraints and state inequality constraints. The method approximates each of the system state variables and each of the control variables by a Bézier curve of unknown control points. The new approximated problems converted to a quadratic programming problem which can be solved more easily than the original problem. Some examples are given to verify the efficiency and reliability of the proposed method.
Mathematical subject classification: 49N10.
Key words: Bézier control points, constrained optimal control, linear time varying dynamical systems.
1 Introduction
In the numerical solution of differential equations, polynomial or piecewisepolynomial functions are often used to represent the approximate solution [1] . Legendre and Chebyshev polynomials are used for solving optimal control problems (see [2] , [3] , [4] and [5] ). Razzaghi and Yousefi [6] defined functions which they called Legendre wavelets for solving constrained optimal control problem. In particular, Bsplines and Bezier curves have become popular tools for solving dynamical systems [7] . For both Bezier and Bspline representations, the curve or surface shape is outlined by a "control polygon" that is formed by connecting the control points. Since the control point structure captures important geometric features of the Bezier or Bspline shape, it is tempting to perform computations just based on the control points. This paper intends to investigate the use of the control points of Bezier representations for solving optimal control problems with inequality constrained. Since the use of Bsplines is to cause the continuity of control curves, we propose to represent the approximate solution in Bezier curves for time varying optimal control problem with inequality constraints. The choice of the Bezier form rather than the Bspline form is due to the fact that the Bezier form is easier to symbolically carry out the operations of multiplication, composition and degree elevation than the Bspline form. We choose the sum of squares of the Bezier control points of the residual to be the measure quantity. Minimizing this quantity gives the approximate solution. Obviously, if the quantity is zero, then the residual function is also zero, which implies the solution is the exact solution. We call this approach the controlpointbased method.
Consider the following time varying optimal control problem with pathwise state inequality constraints,
where
are matrices functions and K(t) = (k_{1}(t) ... k_{p}(t)), R(t) = (r_{1}(t) ... r_{m}(t)), F(t) = (f_{1}(t) ... f_{p}(t))^{T} are vectors functions, where the entries of mentioned matrices and polynomials in [t_{0}, t_{f}] . One may note when the entries of mentioned matrices are not polynomials in [t_{0}, t_{f}] , Taylor series is used for approximation. x(t) = (x_{1}(t) x_{2}(t) ... x_{p}(t))^{T}∈ R^{p} is p ×1 trajectory state vector, u(t) ∈ = {(u_{1}(t) u_{2}(t) ... u_{m}(t))^{T }: < u_{i} < , 1 < i < m} ⊆ R^{m} is m × 1 control vector. The fixed finite terminal time t_{f} is given, and x_{0} is the vector of initial conditions. In addition, l_{i }: R^{p }→ R for i ∈ E, g_{h }: R^{p }→ R for h ∈ I, q : [t_{0}, t_{f}] × R^{p }→ R, and , , 1 < i < m are given bounded functions.
The constrained optimal control of a linear system with a quadratic performance index has been of considerable concern and is well covered in many papers (see [6] and [8] ). One of the methods to solve constrained optimal control problem (1), based on parameterizing the state/control variables, which convert the problem to a finite dimensional optimization problem, i.e. a mathematical programming problem (see [9] , [10] , [11] [16] ).
Analytical techniques developed in [5] are of benefit also in studying the convergence properties of related algorithms for solving optimal control problems, involving Chebyshev type functional constraints where, owing to the use of a variable stepsize in integration or high order integration procedures, it is either not possible or inconvenient to base the analysis on an a priori discretization of the dynamic. The method used slack variables to convert the inequality constraints into equality constraints. This approach has an obvious disadvantage when applied to constrained optimal control of time varying linear problems because it converts the linear constraints into a nonlinearone (see [17] ).
In this paper, we show a novel strategy by using the Bézier curves to find the approximate solution for (1). In this method, we divided the time interval, into k subintervals and approximate the trajectory and control functions in each subinterval by Bézier curves. We have chosen the Bézier curves as piecewise polynomials of degree n, and determine Bézier curves on any subinterval by n + 1 control points. By involving a least square optimization problem, one can find the control points, where the Bézier curves that approximate the action of control and trajectory, can be found as well (see [7] and [18] ).
To show the effectiveness of this method the computational results of threeexamples are presented and compared with the results obtained in [5] and [17] .
2 Least square method
Let k be a chosen positive integer and {t_{0} < t_{1} < ... < t_{k} = t_{f}} be an equidistance partition of [t_{0}, t_{f}] with length τ and S_{j} = [t_{j}_{  1}, t_{j}] for j = 1, 2,..., k. We define the following suboptimal control problems
where C_{jk} = (t_{f})H(t_{f})x_{k}(t_{f}), and δ_{jk} is the Kronecker delta whichhas the value of unity when j = k and otherwise is zero.
are respectively vectors of x(t) and u(t) when are considered in S_{j} = [t_{j}_{  1}, t_{j}] .
Our strategy is using Bezier curves to approximate the solutions x_{j}(t) and u_{j}(t) by v_{j}(t) and w_{j}(t) respectively, where v_{j}(t) and w_{j}(t) are given below. Individual Bezier curves that are defined over the subintervals are joined together to form the Bezier spline curves. For j = 1, 2,..., k, define the Bezier polynomials of degree n that approximate the actions of x_{j}(t) and u_{j}(t) over the interval S_{j} = [t_{j}_{  1}, t_{j}] as follows:
where
is the Bernstein polynomial of degree n over the interval [t_{j}_{  1}, t_{j}] , and are respectively p and m ordered vectors from the control points. By substituting (3) in (2), one may define residual R_{1, j}(t) and R_{2, j}(t) for t ∈ [t_{j}_{  1}, t_{j}] as:
Beside the boundary conditions on v(t), in nodes, there are also continuity constraints imposed on each successive pair of Bezier curves. Since the differential equation is of the first order, the continuity of x (or v) and its first derivativegives
where (t_{j}) is the sth derivative of v_{j}(t) with respect to t at t = t_{j}.
Thus the vector of control points (for r = 0,1, n  1 and n) must satisfy
One may recall that is an p ordered vector. This approach is called the subdivision scheme (or τrefinement in the finite element literature), in the Section 3, we prove the convergence in the approximation with Bezier curves when n tends to infinity.
Note 1: If we consider the C^{1} continuity of w, the following constraints will be added to constraints (5),
where the socalled is an m ordered vector.
Now, we define residual function in S_{j} as follows
where . is the Euclidean norm (Recall that R_{1, j}(t) is a p vector where t ∈ S_{j}) and M is a sufficiently large penalty parameter. Our aim is to solve the following problem over ,
The mathematical programming problem (7) can be solved by many subroutines, we used Maple 12 to solve this optimization problem.
Note 2: If in time varying optimal control problem (1), x(t_{f}) be unknown, then we set C_{kk} = 0.
Note 3: To find good polynomial approximation, one needs to increase the degree of polynomial, in this article, we used sectional approximation, to find accurate results by low degree polynomials.
3 Convergence analysis
In this section, we analyze the convergence of the controlpointbased method when applied to the following time varying optimal control problem:
where the state x(t) ∈ R, u(t) ∈ R, a, b ∈ R, H(t), P(t), Q(t), K(t), R(t), A(t), B(t) and F(t) are polynomials in [0,1] .
Without loss of generality, we consider the interval [0,1] instead of [t_{0}, t_{f}] , since one can change the variable t with the new variable z by t = (t_{f } t_{0}) z + t_{0} where z ∈ [0,1] .
Lemma 3.1. For a polynomial in Bezier form
we have
where a_{i,n}_{1 }_{+ m}_{1} is the Bezier coefficients of x(t) after it is degreeelevated to degree n_{1 }+ m_{1}.
Proof. See [1].
The convergence of the approximate solution could be done in two ways
1. Degree raising case of the Bezier polynomial approximation.
2. Subdivision case of the time interval.
In the following we prove the convergence in each case.
3.1 Degree raising case
Theorem 3.2. If time varying optimal control problem (8) has a unique, C^{1} continuous trajectory solution , C^{0} continuous control solution , then the approximate solution obtained by the controlpointbased method converges to the exact solution as the degree of the approximate solution tends to infinity.
Proof. Given an arbitrary small positive number ε > 0, by the Weierstrass Theorem (see [19] ) one can easily find polynomials Q1,N_{1}(t) of degree N_{1} and Q2,N_{2}(t) of degree N_{2} such that
where ._{∞} stands for the L_{∞}norm over [0, 1] . Especially, we have
In general, Q1,N_{1}(t) and Q2,N_{2}(t) do not satisfy the boundary conditions. After a small perturbation with linear and constant polynomials αt + β and γ, respectively for Q1,N_{1}(t) and Q2,N_{2}(t), we can obtain polynomials P1,N_{1}(t) = Q1,N_{1}(t) + (αt + β) and P2,N_{2}(t) = Q2,N_{2}(t) + γsuch that P1,N_{1}(t) satisfy the boundary conditions P1,N_{1}(0) = a, P1,N_{1}(1) = b, and P2,N_{2}(0) = a_{1}. Thus Q1,N_{1}(0) + β = a, and Q1,N_{1}(1) + α + β = b. By using , one have
Since
so
By the time, a_{1} = P2,N_{2}(0) = Q2,N_{2}(0) + γ, so
Now, we have
Now, let define
for every t ∈ [0, 1] . Thus for N > max{N_{1}, N_{2}}, one may find an upper bound for the following residual:
where C_{1} = 1 + A(t)_{∞} + B(t)_{∞} is a constant.
Since the residual R(P_{N}) : = LP_{N}(t)  F(t) is a polynomial, we can represent it by a Bezier form. Thus we have
Then from Lemma 1 in [1] , there exists an integer M(> N) such that when m_{1} > M, we have
which gives
Suppose x(t) and u(t) are approximated solution of obtained by the controlpointbased method of degree m_{2} (m_{2}> m_{1}> M). Let
Define the following norm for difference approximated solution (x(t), u(t)) and exact solution ((t), (t)):
It is easy to show that:
Last inequality in (12) is obtained from Lemma 1 in [1] in which C is a constant positive number. Now from Lemma 1 in [1] , one can easily show that:
where last inequality in (13) is coming from (10).
Thus, from (13) we have:
Since the infinite norm and the norm defined in (11) are equivalent, there is a ρ_{1} > 0 where
Now, we show that the approximated cost function tends to exact cost function as the degree of Bezier approximation increases. Define
for t ∈ [0, 1] . Now, there are four positive integers M_{i} > 0, i = 1,..., 6, such that
Since
we have
so
now, we have
This completes the proof.
3.2 Subdivision case
Theorem 3.3 Let (x, u) be the approximate solution of the linear optimal control problem obtained by the subdivision scheme of the controlpointbased method. If has a unique solution (, ) where (, ) is smooth enough so that the cubic spline T(, ) interpolating to (, ) converges to (, ) in the order O(τ^{q}), (q > 2), where τ is the maximal width of all subintervals, then (x, u) converges to (, ) as τ → 0.
Proof. We first impose a uniform partition Π_{d} = ∪_{i}[t_{i}, t_{i}_{ + 1}] on the interval [0, 1] as t_{i} = id where d =
Let be the cubic spline over Π_{d} interpolating to (, ). Then for an arbitrary small positive number ε > 0, there exists a δ_{1} > 0 such that
provided that d < δ_{1}. Let
be the residual. For each subinterval [t_{i}, t_{i}_{ + 1}] , is a polynomial. On each interval [t_{i}, t_{i}_{ + 1}] , we impose another uniform partition Π_{i, }τ = ∪_{j} [t_{i, j}, t_{i , j }_{+ 1}] as t_{i,j} = id + jτ where τ = , j = 0,..., m_{1}. Express in [t_{i, j }_{ 1}, t_{i, j} ] as
By Lemma 3 in [1] , there exists a δ_{2} > 0 (δ_{2}< δ_{1}) such that when τ < δ_{2}, we have
Thus
or
Now combining the partitions Π_{d} and all Π_{i,}_{τ} gives a denser partition with the length τ for each subinterval. Suppose (x(t), u(t)) is the approximate solution by the controlpointbased method with respect to this partition, and denote the residual over [t_{i,j}_{  1}, t_{i, j}] by
Then there is a constant C such that
last inequality in (14) is obtained from Lemma 1 in [1] . One can easily show that:
Thus, from (15) we have:
Since the infinite norm and the norm defined in (11) are equivalent, there is a ρ_{2} > 0 where
Now, we show that the approximated cost function tends to exact cost function as the degree of approximation increases. Define
for t ∈ [t_{i,j}_{  1}, t_{i,j}] . Now, there are four positive integers M_{i} > 0, i = 1,..., 6, such that P(t)_{∞} < M_{1}, Q(t)_{∞} < M_{2}, K(t)_{∞} < M_{3}, R(t)_{∞} < M_{4}, (t)_{∞} < M_{5}, and (t)_{∞} < M_{6}. Since
we have
so
now, we have
Now, from Lemma 3 in [1] , we conclude that the approximated solution converges to the exact solution in order o(τ^{q}), (q > 2).
This completes the proof.
4 Numerical examples
Example 4.1. Consider the following constrained optimal control problem (see [5] ):
Let k = 6 and, n = 3, if one consider the C^{1} continuity of u(t), one can find the following solution from (7)
Then by involving the above control function u(.), in the indicated differential equation, one can find the trajectories x_{1}(.) and x_{2}(.) as follows:
and,
The graphs of approximated trajectories x_{1}(t) and x_{2}(t) are shown respectively in Figure 1 and Figure 2, and the graph of approximated control is shown Figure 3. The approximated and exact objective function are respectively I = 5.389857789, I^{*} = 5.528595476 (see [5] ).
In this example, we used Bezier polynomials of degree 10 to approximate the trajectories x_{1}(.) and x_{2}(.) through the time interval J = [0,3] , and without using subintervals. The objective function is found I = 5.360252730 and the trajectories x_{1}(t), and x_{2}(t) for t ∈ J are shown in Figures 4, 5. For reduction in complicated manipulations, the authors recommend to use low degree Bezier approximation polynomials and use subintervals approach.
Example 4.2. Consider the following optimal control problem ([17] ):
Let k = 12, t_{0} = 0, t_{f} = 1, t_{j} = t_{0}+ (j = 1,..., 12), and n = 3. From , one can find the following solution
The graphs of approximated trajectories x_{1}(t) and x_{2}(t) are shown respectively in Figure 6 and Figure 7, and the graph of approximated control is shown Figure 8. The approximated and exact objective function respectively are I = 0.1728904530, I^{*} = 0.17078488 (see [17] ). The computation takes 5 seconds of CPU time when it is performed using Maple 12 on an AMD Athelon X4 PC with 2 GB of RAM. The QPSolve command solves (7), which involves computing the minimum of a quadratic objective function possibly subject to linear constraints. The QPSolve command uses an iterative active set method implemented in a builtin library provided by the numerical algorithms group.
5 Conclusions
In this paper, a Bézier control points method for solving optimal control problems governed by time varying linear dynamical systems with terminal state constraints and inequality constraints on the states and control has been suggested. The method replaces the constrained optimal control problem by a quadraticprogramming one, and the control point structure provides a bound on the residual function. The new mathematical programming problem is intuitive and easy to solve.
Numerical examples show that the proposed method is reliable and efficient. The control polygon gives an intermediate approximation to the appearance of the solution.
Acknowledgment. The authors like to express their sincere gratitude to referees for their very contractive advises.
Received: 01/VIII/11.
Accepted: 13/X/12.
#CAM396/11.
 [1] J. Zheng, T.W. Sederberg and R.W. Johnson, Least squares methods for solving differential equations using Bezier control points Appl. Num. Math, 48 (2004), 237252.
 [2] M. Razaghi and G. Elnagar, A legendre technique for solving timevarying linear quadratic optimal control problems J. Franklin Institute, 330 (1993), 453463.
 [3] G.N. Elnagar and M. Razzaghi, A Collocationtype method for linear quadratic optimal control problems Optimal. Control. Appl. Methods, 18 (1997), 227235.
 [4] L.A. Rodriguez and A. Sideris, A sequential linear quadratic approach for constrained nonlinear optimal control American Control Conference, San Francisco, CA, USA (2011).
 [5] J. Vlassenbroeck, A Chebyshev polynomial method for optimal control with state constraints Automatica, 24 (1988), 499504.
 [6] M. Razaghi and S. Yousefi, Legendre wavelets method for constrained optimal control problems Mathematics Methods in the Applied Sciences, 25 (2002), 529539.
 [7] M. Gachpazan, Solving of time varying quadratic optimal control problems by using Bézier control points Computational & Applied Mathematics, 30(2) (2011), 367379.
 [8] V. Yen and M. Nagurka, Optimal control of linearly constrained linear systems via state parameterization Optimal Control. Appl. Methods, 13 (1992), 155167.
 [9] G. Elnagar, M. Kazemi and M. Razzaghi, The pseudospectral Legendre methodfor discretizing optimal control problem IEEE Trans. Automat. Control., 40 (1965), 17931796.
 [10] P.A. Frick and D.J. Stech, Solution of optimal control problems on a parallel machine using the Epsilon method Optimal Control. Appl, 16 (1995), 117.
 [11] C.P. Neuman and A. Sen, A suboptimal control algorithm for constrained problems using cubic splines Automatica, 9 (1973), 601613.
 [12] R. Pytlak, Numerical methods for optimal control problems with state constraints First ed, Berlin: SpringerVeriag (1999).
 [13] H.R. Sirisena, Computation of optimal controls using a piecewise polynomial parameterization IEEE Trans. Automat. Control, 18 (1973), 409411.
 [14] H.R. Sirisena and F.S. Chou, State parameterization approach to the solution of optimal control problems Optimal Control. Appl. Methods, 2 (1981), 289298.
 [15] K. Teo, C. Goh and K. Wong, A unified computational approach to optimal control problems First edition, London: Longman, Harlow (1981).
 [16] I. Troch, F. Breitenecker and M. Graeff, Computing optimal controls for systems with state and control constraints, in: Proceedings of the IFAC Control Applications of Nonlinear Programming and Optimization, (1989), 3944.
 [17] H. Jaddu, Spectral method for constrained linearquadratic optimal controlMath. Comput. in Simulation, 58 (2002), 159169.
 [18] M. Evrenosoglu and S. Somali, Least squares methods for solving singularity perturbed twopoints boundary value problems using Bézier control point Appl. Math. Letters, 21(10) (2008), 10291032.
 [19] W. Rudin, Principles of mathematical analysis, McGrawHill (1986).
Publication Dates

Publication in this collection
28 Nov 2012 
Date of issue
2012
History

Received
01 Aug 2011 
Accepted
13 Oct 2012