OPTIMAL PROGRAM CONTROL IN THE CLASS OF QUADRATIC SPLINES FOR LINEAR SYSTEMS Текст научной статьи по специальности «Математика»

2020 ВЕСТНИК САНКТ-ПЕТЕРБУРГСКОГО УНИВЕРСИТЕТА Т. 16. Вып. 4

ПРИКЛАДНАЯ МАТЕМАТИКА. ИНФОРМАТИКА. ПРОЦЕССЫ УПРАВЛЕНИЯ

ПРОЦЕССЫ УПРАВЛЕНИЯ

UDC 517.977.5 MSC 49N10

Optimal program control in the class of quadratic splines for linear systems*

A. S. Popkov

St. Petersburg State University, 7—9, Universitetskaya nab., St. Petersburg, 199034, Russian Federation

For citation: Popkov A. S. Optimal program control in the class of quadratic splines for linear systems. Vestnik of Saint Petersburg University. Applied Mathematics. Computer Science. Control Processes, 2020, vol. 16, iss. 4, pp. 462-470. https://doi.org/10.21638/11701/spbu10.2020.411

This article describes an algorithm for solving the optimal control problem in the case when the considered process is described by a linear system of ordinary differential equations. The initial and final states of the system are fixed and straight two-sided constraints for the control functions are defined. The purpose of optimization is to minimize the quadratic functional of control variables. The control is selected in the class of quadratic splines. There is some evolution of the method when control is selected in the class of piecewise constant functions. Conveniently, due to the addition/removal of constraints in knots, the control function can be piecewise continuous, continuous, or continuously differentiable. The solution algorithm consists in reducing the control problem to a convex mixed-integer quadratically-constrained programming problem, which could be solved by using well-known optimization methods that utilize special software.

Keywords: optimal control, differential equations, linear control system, quadratic spline, mixed-integer quadratically-constrained programming.

1. Introduction. Often, in mathematical models describing various mechanical, electrodynamic, economic, or other processes, the control function by its nature must be continuous or even smooth in time. Another example is when smooth movement is an important goal of the control problem. Thus, in [1] the problem of control constructing for a wheeled mobile robot is considered. Such machines are used to transport people with disabilities, so it is necessary to do this with comfort for passengers.

In such situations, it is convenient to select the control from the class of quadratic splines, because these functions are quite easy to research, and they can approximate difficult nonlinear functions. Also by adding or removing constraints in knots, the control

* The reported study was funded by the Russian Foundation for Basic Research (project N 19-3190033).

© Санкт-Петербургский государственный университет, 2020

function can be made piecewise-continuous, continuous, or smooth. The entire considered time segment is splitted into some predetermined number of sub-segments N (for simplicity, we assume that all sub-segments are of the same length, but this is easy to change). And on each segment, we will consider a polynomial of the second degree as a control function.

This work is a continuation of the articles [2,3], in which algorithms of constructing optimal control for the linear and nonlinear cases in the class of piecewise constant functions are described, which is very justified, given the global digitalization and, therefore, the sampling of most processes. The original idea in this area belongs to the R. Gabasov research team [4]. Related researches are also presented in [5,6].

2. Problem statement. The dynamics of the process is consistent with a linear system of ordinary differential equations with control:

x = Ax + Bu, (1)

where x(t) is a n-dimensional phase vector; X(t) is a derivative of x(t); u(t) is a r-di-mensional control vector; A is a (n x n)-matrix; B is a (n x m)-matrix; t is a time variable defined in a segment [0, T].

The terminal problem is to transfer the system from the initial to the final state:

x(0) = x*, x(T) = x*. (2)

The controls have two-sided constraints:

< Ui < lf\ ¿ = 1~F. (3)

We also assume that the controls belong to the class of smooth functions:

Ui(t) £ C1 ([0, T]), (4)

The objective of the control is to minimize the integral of the quadratic form of u:

T

J = J u(t)TQu(t) dt —> min . (5)

0

In this expression, Q is a m x m symmetric matrix. Suppose that Q is a Positive Semi-Definite matrix. This is an optional restriction. But the convexity of the problem depends on this, which is very important for solvability.

The control is chosen in the class of uniform quadratic splines with knots tk = Tk/N, k = 0, N, where N is predefined number:

, Pikt2+Pikt+Pil\ iefa-i, tk), .

i{t) = \ __i=l,r. (6)

k = l,N,

Thus the problem is in finding optimal admissible values of pf^, Pik , Pik .

3. Reduction to MIQCP problem. In this section it will be shown that the optimal control problem (1)-(6) can be reduced to an optimization problem with continuous and binary variables, linear and quadratic constraints and a quadratic objective function. This problem termed a Mixed-Integer Quadratically-Constrained Programming (MIQCP).

Smoothness condition. Since the control is represented as a quadratic spline, there are guaranteed smoothness between knots. Therefore, it remains to ensure continuity and smoothness at knots.

Continuity at knots for control functions:

Pik4+p[S+Pik=P^+A+P^+itk+p[Zv k=T^T, * = (7)

Continuity at knots for derivatives of control functions:

k = i = (8)

Note, if smoothness is not required at knots, we just need to disable restrictions (8). If continuous is also not required, we also need to disable restrictions (7).

Straight constraints. In accordance with inequalities (3) and representation of the control functions as (6), the following conditions must be satisfied:

^ < P$t2 + < f'1 Vi € [ifc-1, tk), k = t =

These conditions can be transformed to extremum conditions

l(1} < tp mf 2 + + Pk) ,

tE[tk-1, tk ) V /

^ > sup (p^t2 + pkt + pk) ,

te[tk-i, tk) v 7

k=l,N, i = l,r.

In this case, on each segment [tk~i, tk), it is necessary to consider only three points suspicious of extremum: tk-1, tk and —p^/p^. There are simple restrictions for knot points:

¿(1) < p(2)t2 + p(l)t. 1 + p(0) < ,(2) li < pik bk-1T pik °k-1T pik < ,

,(1) < p(2)t2 + p(1)t + p(0) < ¿(2) li < pik tk + pik k + pik < li ,

k=l,N, i = l,r. (9)

More difficult situation with the parabola extremum — p^kk /2pfk. It would be wrong to say that the extremum value must satisfy two-sided boundaries, because, after all, it is permissible to violate the boundaries if the argument of the extremum is to the left or right of the segment [tk-i, tk).

For simplicity reassign some variables:

P := pik, Ti := tk-i, T2 := tk, li := l^, ¡2 := ¡f^ And suppose that values of quadratic polynomial at knots are fixed:

y := P(2)t2 + p(i)t + p(0) y := p(2)t2 + p(i)t + p(0) Vi := Pik tk-i + Pik tk-i + Pik , V2 := Pik tk + Pik tk + Pik ■

Then it is possible to express the values of p^k and pf° through p, yi and y2:

(i) V2 — yi , S (0) V2 — Vi Pik =--P(T 1 + T2j, Pik =PTlT2--n + y 1.

T2 — Ti T2 — Ti

Now we can find the conditions that an extremum value of quadratic polynomial is equal l2 :

V2 -yi , , A2 . ( V2 - VI , n

- P(ti + r2 j -4p\ priT-2--ri +2/1 — ¿2 =0.

T - T1 ) \ T2 - T1

After transformations we get the quadratic equation of p:

(i2 - il) V + (—2(2/1 + y2) + 4/2)p + fyi^Jtl) = 0,

V T2 - T1 )

above is a quadratic equation for unknown variable p that has two solutions:

f Vh -yi± Vh ~ V2 P = - -

V T2 - T1

Thus, the task was solved: find the parameters of the quadratic polynomial, when two points are fixed, and the extremum is equal to a predeterminated value. Analysis of solutions makes it easy to understand that the solution with a larger absolute value corresponds to the extremum value located between two fixed points, and another solution corresponds to the extremum value located outside these points.

By doing the same for the lower bound, we get a two-sided constraint for p:

_ _ 2 _ _ 2

(Vh -y 1 + aA — 2/2A <p< f Vyi ~ll + -/A ^^

V T~2 -n J ^ ^ V T2 - ri J '

Now let's return to the previous notations and introduce new variables s^, s(2\ s(j), s^k and add restrictions for them:

(1) _ (2),2 + p(l), + p(0) l(l)

sik _ Pik tk-1 + Pik tk-1 + Pik - li ,

s(2) _ p(2)t2 + p(1)t + p(0) _ l(1) sik pik bk> pik bk "I" pik li ,

(3) _ /(2) _ (2) 2 _ (1) (0) k = l^N, i=I7F. (11)

sik li pik bk-1 pik bk-1 pik , ' ' ' v '

s(4) _ l(2) _ p(2)t2 _ p(1)t p(0) sik li pik Lk pik tk pik ,

sk £ 0, sk £ 0, sk £ 0, sk £ 0, ik ik ik ik

Then rewrite equation (10):

k=l,N, i = l,r. (12)

Next step is to split variables p^2 to new non-negative variables r^, r(2, r^, r^:

t 2P2 = N 2 ( - ri -Г2++rj?)

k=l,N, i = l,r. (13)

r(1) > о r(2) > о r(3) > о r(4) > о r ik > r ik > r ik > r ik >

The idea behind this separation is that in equations (12) r(1 should be responsible for the

(2)

linear part of the lower boundary, r\k' should be responsible for the non-linear part of the

(3)

lower boundary, r\k should be responsible for the linear part of the upper boundary and

r(k should be responsible for the non-linear part of the upper boundary.

(2) (1) (2) (2) Important detail is that if p\k' > 0 then r\k + r\k' = 0. Otherwise, if p\k' < 0 then

r(3) + rk = 0. To ensure this, we have to introduce a binary variable Sik and append following constraints:

v J k=l,N, i = l,r. (14)

ri3) + r|4) < 4 (f> - lf>) (1 - Sik)

Thereby, if Sik = 0 then rH*1 + rij) = 0, else if Sik = 1 then rf3 + r^ = 0. Element 4(l¿2^ — l¿^) means upper bound because it is easy to demonstrate that if pi2 > 4(l|2^ — l¿^)N2/T2 then straight constraints are anyway violated.

Finally, let's rewrite the linear part of constraints (12) in the new notations:

W^sW + stf, rfk^s^ + s^, k = Tjf, ¿ =T"F, (15)

and non-linear part:

Note that constraints (16) are a rotated Second-Order Cone constraints. This is main achievement and we will discuss it in the Section 4.

In summary, it was shown how two-sided constraints (3) were transformed into linear constraints (9), (11), (13), (15), linear inequalities with binary variables (14) and quadratic inequalities (16).

Terminal conditions. The condition for the initial state is already satisfied, since it is a part of the Cauchy problem. Let us define the conditions for the final state. Write the Cauchy formula at final point and apply conditions (2)

T

eATx* +У eA(T-t)Bu(t)dt.

Let's divide segment [0, T] into N subsegments with the same lengths and apply that the control function is selected in the form (6)

e x^ -+

r N

+EE

i=1 k=1

( tk \ ( tk \ ( tk

J e^^dt^pk + J eA(T-t)bltdt^pi1k> + J eA(T-t)bidt^p^} \tk-1

e

i i \tfc-i i \tfc-i

\tfc-i

Here bi denotes the column i of the matrix B.

*

x

:

x

Vector functions eA(T t)bi could be numerically calculated as solutions of the following Cauchy problems (it's profitable if r < n):

¿« = -AzW(t), . _

i = 1, r.

z(i)(T ) —

Now let's introduce new notations:

__* AT

g .— x e x^,

di'k'0 :— f eA(T-t)bidt, ifc-i

d^'1 := / eA(T~%itdt, k = h~N, i = 1~F, tk-l

di'k'2 :— f eA(T-t)bit2dt, tk-i

it is assumed that integrals will be calculated numerically.

Then we get n linear equalities

N

,k,0(0). >,M (1) , Ai,k,2(2)^

E E + +VMV) •■/;• (it)

i=1 k=1

Objective function. Replace the control u(t) by quadratic spline (6) in the objective function (5):

N r r

J = NTT j {p(2k t2 + p(k t+pty t2+pCk t+pSI) dt,

k=1 il =1 i2 = 1t

t k-l

where qili2 is an element of the matrix Q.

After calculating of the integrals and introducing new designations, we obtain quadratic form in the variables pf°, p(k and pf^:

N r r 2 2

j — EEEE Ehia^pao, (18)

k=1 i1 = 1 i2 = 1 a1=0 a2 = 0

here

1 T P / \

KXak ■■= r - - ln> /3 = «1 + «2 +1,

k=l,N, i\ = 1, r, ¿2 = l,r, ai=0,2, «2=0,2.

According Q is a positive semi-definite matrix, in 3 steps it can be proved that the quadratic form in the formula (18) is a positive semi-definite:

1) the resulting matrix is a block diagonal matrix consisting of N blocks for each к = 1 ,N. Therefore, it is necessary to demonstrate that any block is a positive semi-definite matrix;

2) for each k if we fix t then ui(t) is a linear function of p(0), p(k and pfk • Next it can be shown that composition G(F(£)) is a positive semi-definite quadratic form, if G(n) is a positive semi-definite quadratic form and F(£) is a linear function. Thence we see that under the integral is a positive semi-definite quadratic form in the variables p(0), pk and Pk for each fixed t;

3) the last step is to show that the result of integrating of a positive semi-definite quadratic form with the parameter t is a positive semi-definite quadratic form.

4. Model analysis. In previous section the optimal control problem (1)-(6) was transformed to the optimization problem (7)-(9), (11), (13)-(18). In this article we are not aimed in receiving matrix-based model with a minimal number of variables. Because it is quickly achieved at the stage of presolve in the runtime of the optimization software. Real purpose is to obtain effective model, optimal in the terms of calculation time. For ease of narration, we have described the model in a constraint-based form. So let's take a closer look at this.

Variables:

1) 3rN basic continuous variables: pfk , p^, pfk , k = 1, N, i = 1, r.

2) 8rN continuous non-negative variables: s\l\ sfk , sfk , s^, rfj), k = I, N, i = T~r. _

3) rN binary variables: Sik, k = 1, N, i = 1, r.

Constraints:

1) 7rN — 2r + n linear equalities: (7), (8), (11), (13), (17). So a number of variables can be excluded.

2) 6rN linear inequalities: (9), (15).

3) 2rN linear inequalities with binary variables: (14).

4) 2rN quadratic inequalities: (16). Although quadratic constraints are not represented by positive semi-definite quadratic forms, they are rotated Second-Order Cone equations. And it can be proved that in the case when the variables .s(1>, sfk , s^, s^ are non-negative, there is a convex set of solutions.

Accordingly, we have a convex set of solutions in the optimization problem obtained by linear and quadratic constraints.

Objective. The objective function (18) of the optimization problem is a positive semi-definite quadratic form in the variables pfk , p^, pfk , k = 1 ,N, i = 1, r.

Summary. Thereby, we have obtained the optimization problem with continuous and binary variables, linear and convex quadratic constraints, and a convex quadratic objective. According to the generally accepted classification, this problem is a convex case of MIQCP. There are various software that can resolve this problem, such as Gurobi, CPLEX, Artelys Knitro, etc. [7].

5. Example. The problem of the damping of the sleeping Lagrange top [8] was considered as an illustration of the method:

| + Bn — = ui, r/ — B£ — An = u2,

here £(t) and n(t) are denote coordinates of the top end; ui(t), u2(t) are generalized forces; A and B are model parameters, indicating the physical properties of the object. Let's put A =1, B = 3.

The terminal problem is to change object position (£, n) from ( — 10,5, 6,8) to (0,0,0,0). The value of T was set to 20, N also equals 20. All lower bounds are —2, all upper bounds are 2. And finally, Q is two-dimensional identity matrix.

MIQCP problem was obtained and then solved using Gurobi. Total processing time is 1.56 s. The result is illustrated at Figure.

иi(t), uAt) Control Function

rj(t) Phase Portrait

(moving an object from point (—10, 6) to point (0, 0). The object's speed also fades to zero)

6. Conclusion. Recall that the initial problem is a control problem with a linear ordinary differential equations (ODE) system, a quadratic objective functional, fixed initial and final states and bounded control in the form of quadratic splines. This problem is a some alternative or generalization of similar problem with a control in the form of piece-wise-constant function. The optimal control problem has been reduced to the convex case of MIQCP. In order to accomplish this, it is necessary to carry out a numerical calculation of ODEs and definite integrals. The main achievement of current research that the derived optimization problem is a convex, which means that the global extremum can be efficiently calculated using special software.

This study is completed, but further researches may be related to the consideration of nonlinear ODE systems or consideration of other control functions instead of quadratic splines. More specific restrictions to the phase or control variables may also be considered.

References

1. Park J. J., Kuipers B. A smooth control law for graceful motion of differential wheeled mobile robots in 2D environment. 2011 IEEE International Conference on Robotics and Automation, 2011, pp. 4896-4902.

2. Popkov A. S. Application of the adaptive method for optimal stabilization of a nonlinear object. 2016 International Conference Stability and Oscillations of Nonlinear Control Systems (Pyatnitskiy's Conference), 2016, pp. 1-3.

3. Popkov A. S., Smirnov N. V., Smirnova T. E. On modification of the positional optimization method for a class of nonlinear systems. ACM International Conference Proceeding Series, 2018, pp. 46-51.

4. Balashevich N.V., Gabasov R., Kirillova F. M. Chislennye metody programmnoy i pozitsionnoy optimizatsii lineynyh sistem upravleniya [Numerical methods of program and positional optimization

of linear control systems]. Journal of Computational Mathematics and Mathematical Physics, 2000, vol. 40(6), pp. 838-859. (In Russian)

5. Girdyuk D.V., Smirnov N. V., Smirnova T. E. Optimal control of the profit tax rate based on the nonlinear dynamic input-output model. ACM International Conference Proceeding Series, 2018, pp. 80-84.

6. Baranov O.V., Smirnov N.V., Smirnova T. E., Zholobov Y. V. Design of a quadrocopter with PID-controlled fail-safe algorithm. Journal of Wireless Mobile Networks, Ubiquitous Computing, and Dependable Applications, 2020, vol. 11(2), pp. 23-33.

7. GAMS Development Corporation: Solver Manuals. Available at: http://www.gams.com/latest/docs/S_MAIN.html#SOLVERS_MODEL_TYPES (accessed: September 15, 2020).

8. Babadzanjanz L. K., Pototskaya I. Yu. Upravlenie po kriteriyu raskhoda v mehanicheskih sistemah [Control by expense criteria in mechanical systems]. St. Petersburg, Saint Petersburg State University Press, 2003, 137 p. (In Russian)

Received: October 08, 2020. Accepted: October 23, 2020.

Authors' information:

Alexander S. Popkov — Postgraduate Student, Research Engineer; alexandr.popkoff@gmail.com

Оптимальное программное управление в классе квадратичных сплайнов для линейных систем*

А. С. Попков

Санкт-Петербургский государственный университет, Российская Федерация, 199034, Санкт-Петербург, Университетская наб., 7-9

Для цитирования: Popkov A. S. Optimal program control in the class of quadratic splines for linear systems // Вестник Санкт-Петербургского университета. Прикладная математика. Информатика. Процессы управления. 2020. Т. 16. Вып. 4. С. 462-470. https://doi.org/10.21638/11701/spbu10.2020.411

Разработан алгоритм решения задачи оптимального управления в случае, когда рассматриваемый процесс описывается линейной системой обыкновенных дифференциальных уравнений, заданы начальное и конечное состояния управляемого объекта и присутствуют двусторонние ограничения на управления. Целевой функционал задан в виде интеграла от положительно полуопределенной квадратичной формы от управлений. Управление выбирается в классе квадратичных сплайнов. При таком виде управлений оказывается весьма удобно, что, благодаря добавлению или удалению ограничений в узлах, функция управления может быть кусочно-непрерывной, непрерывной или гладкой. Алгоритм решения заключается в сведении задачи управления к задаче выпуклого смешанного целочисленного квадратичного программирования с квадратичными ограничениями, которую можно решить при помощи известных методов оптимизации с использованием специального программного обеспечения.

Ключевые слова: оптимальное управление, дифференциальные уравнения, линейная управляемая система, квадратичные сплайны, смешанное целочисленное квадратичное программирование с квадратичными ограничениями.

Контактная информация:

Попков Александр Сергеевич — аспирант, инженер-исследователь; alexandr.popkoff@gmail.com

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (проект № 19-31-90033).