* ==> opt_cont.gms
$title An Optimal Control ELQP
$offsymxref offsymlist offuellist offuelxref
$inlinecom /* */
/* %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
This file contains an Extended Linear-Quadratic Program (ELQP)
formulated as an MCP. See
R.T. Rockafellar, "Linear-Quadratic Programming and Optimal
Control", SIAM J. Control and Optimization 25 3, 1987, p. 781-814.
In general, an ELQP consists of primal and dual problems of the
following form:
(P): min pu + 1/2 uPu + rho_{V,Q}(q - Ru)
s.t. u in U
(D): max qv - 1/2 vQv - rho_{U,P}(R^v - p)
s.t. v in V
where P and Q are positive semidefinite matrices,
p and q are vectors, u and v are variables in the
feasible sets U \subset R^n, V \subset R^m,
and rho_{X,D} is a piecewise linear-quadratic penalty
function depending on the feasible set X and matrix D.
The optimality conditions for the ELQP are equivalent to
VI(T,U \cross V), where T(u,v) follows. It is this VI which
we formulate as a MCP.
T(u,v) := (Pu - R'v + p)
(Ru + Qv - q)
The ELQP in question arises from the discretization of the
following continuous (primal) optimal control problem:
t_R
integral [p(t)*u(t) + 1/2 u(t)*P(t)*u(t) - c(t)*x(t)] dt +
t_L
p^L*u^L + 1/2 u^L*P_L*u^L - c^R*x(t_L) +
min
t_R
integral rho(q(t) - C(t)*x(t) - D(t)*u(t)) dt +
t_L
rho(q^R - C_R*x(t_R))
s.t. (dx/dt)(t) = A(t)*x(t) + B(t)*u(t) + b(t),
x(t_L) = B_L*u^L + b^L
where u^L are some initial control variables,
u(t) the controls at time t, x(t) the state at time t.
The dual, not given here, has state variables y(t),
controls v(t), and a terminal control variable v_T.
For any N >= 1, the interval [L, R] is divided into N pieces by a
uniform mesh of N+1 points, including endpoints L and R. For the
primal problem, our discretization includes N+1 state vectors:
x^1=x(L), .., x^N, xR = x(R). The control vectors at the first N mesh
points are u^1, .., u^N. No control vector is needed for the right
endpoint; however, there is the control vector u^L==uI,
the initial control which determines state x(t_L).
It is not necessary to specify the dual problem; it is implicit in
the primal. FYI, our dual state variables are yL=y(L), y^1, .., y^N =
y(R), while there is no dual control for the left endpoint; instead,
v^1, .., y^N are the dual control variables at the last N mesh points.
There is an additional control vector yT to enforce a terminal
condition on the dual states.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */
sets
mesh / 1 * 8 /, /* points 1 .. N in primal & dual mesh */
controls / 1 * 8 /, /* controls */
icontrols / 1 * 8 /, /* initial controls */
states / 1 * 8 /, /* states */
/* ds duals to states */
duals / 1 * 8 /, /* dual controls, or constraint multipliers */
tduals / 1 * 8 /; /* terminal dual controls */
option seed = 87;
alias (ds, states),
(controls, controls2),
(icontrols, icontrols2),
(duals, duals2),
(tduals, tduals2);
scalar
pcount, /* count of primal variables */
dcount, /* count of dual variables */
timespan, /* t_R - t_L */
N; /* number of grid (mesh) points */
timespan = 1;
N = card(mesh);
pcount = N * (card(controls) + card(states))
+ card(icontrols) +card(states);
dcount = N * (card(duals) + card(ds))
+ card(tduals) +card(ds);
/* problem size for final MCP is
N * (#controls + #duals + 2 * #states))
+ #icontrols + #tduals + 2 * #states */
parameters
u_lo(controls),
u_up(controls),
uI_lo(icontrols),
uI_up(icontrols),
v_lo(duals),
v_up(duals),
vT_lo(tduals),
vT_up(tduals),
P(controls,controls2),
pp(controls),
PI(icontrols, icontrols2),
ppI(icontrols),
Q(duals,duals2),
qq(duals),
QT(tduals, tduals2),
qqT(tduals),
A(ds, states),
B(ds, controls),
bb(ds),
BI(ds, icontrols),
bbI(ds),
C(duals, states),
cc(states),
CT(tduals, states),
ccT(states),
D(duals, controls);
/* %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
All the parameters have been declared; we must now define
them.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */
u_up(controls) = uniform(-5, 5);
u_lo(controls) = u_up(controls) - uniform(0, 4.5);
uI_up(icontrols) = uniform(-5, 5);
uI_lo(icontrols) = uI_up(icontrols) - uniform(0, 4);
v_up(duals) = uniform(-5, 5);
v_lo(duals) = v_up(duals) - uniform(0, 4);
vT_up(tduals) = uniform(-5, 5);
vT_lo(tduals) = vT_up(tduals) - uniform(0, 4);
P(controls,controls) = uniform(1, 3);
pp(controls) = uniform(-2, 2);
PI(icontrols, icontrols) = uniform(1, 3);
ppI(icontrols) = uniform(-2, 2);
Q(duals,duals) = uniform(1, 3);
qq(duals) = uniform(-2, 2);
QT(tduals, tduals) = uniform(1, 3);
qqT(tduals) = uniform(-2, 2);
A(ds, states) = uniform(-0.9, 0.9);
B(ds, controls) = uniform(-0.2, 0.2);
bb(ds) = uniform(-0.5, 0.5);
BI(ds, icontrols) = uniform(-0.2, 0.2);
bbI(ds) = uniform(-0.5, 0.5);
C(duals, states) = uniform(-0.2, 0.2);
cc(states) = uniform(-0.5, 0.5);
CT(tduals, states) = uniform(-0.2, 0.2);
ccT(states) = uniform(-0.5, 0.5);
D(duals, controls) = uniform(-10, 10);
/* %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
The parameters are defined. Now
declare the variables and equations,
and define the functions and the
box constraints.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */
variables
uI(icontrols), /* initial control variables */
u(mesh, controls), /* control variables for [L,R) */
x(mesh, states), /* state variables for [L,R) */
xR(states), /* state variables at t=R */
v(mesh, duals), /* dual control vars, (L,R] */
vT(tduals), /* terminal dual control vars */
yL(ds), /* dual state vars at t=L */
y(mesh, ds); /* dual state vars, (L,R] */
uI.lo(icontrols) = uI_lo(icontrols);
uI.up(icontrols) = uI_up(icontrols);
u.lo(mesh,controls) = u_lo(controls);
u.up(mesh,controls) = u_up(controls);
v.lo(mesh,duals) = v_lo(duals);
v.up(mesh,duals) = v_up(duals);
vT.lo(tduals) = vT_lo(tduals);
vT.up(tduals) = vT_up(tduals);
/*
uI.lo(icontrols) = (uI_lo(icontrols) + uI_up(icontrols)) / 2;
u.l(mesh,controls) = (u_lo(controls) + u_up(controls)) / 2;
v.lo(mesh,duals) = (v_lo(duals) + v_up(duals)) / 2;
vT.lo(tduals) = (vT_lo(tduals) + vT_up(tduals)) / 2;
*/
equations
deluI(icontrols),
delu(mesh, controls),
delx(mesh, states),
delxR(states),
delv(mesh, duals),
delvT(tduals),
delyL(ds),
dely(mesh, ds);
deluI(icontrols) ..
sum(icontrols2, PI(icontrols, icontrols2) * uI(icontrols2))
+ ppI(icontrols)
- sum(ds, yL(ds) * BI(ds, icontrols))
=g= 0;
delu(mesh, controls) ..
( sum(controls2, P(controls,controls2) * u(mesh,controls2))
+ pp(controls)
- sum(duals, v(mesh,duals) * D(duals,controls))
- sum(ds, y(mesh,ds) * B(ds, controls)) * timespan
) / N
=g= 0;
delx(mesh, states) ..
y(mesh-1,states) + (yL(states)) $(ord(mesh) eq 1)
=e=
y(mesh,states)
+ ( cc(states)
+ sum(ds, y(mesh,ds) * A(ds,states)) * timespan
+ sum(duals, v(mesh,duals) * C(duals, states))
) / N;
delxR(states) ..
- ccT(states)
- sum(tduals, vT(tduals) * CT(tduals, states))
+ sum(mesh$(ord(mesh) eq N), y(mesh,states))
=e= 0;
delv(mesh,duals) ..
( sum(controls, D(duals,controls) * u(mesh,controls))
+ sum(states, C(duals,states) * x(mesh,states))
+ sum(duals2, Q(duals, duals2) * v(mesh,duals2))
- qq(duals)
) / N
=g= 0;
delvT(tduals) ..
sum(states, CT(tduals,states) * xR(states))
+ sum(tduals2, QT(tduals,tduals2) * vT(tduals2))
- qqT(tduals)
=g= 0;
delyL(ds) ..
sum(icontrols, BI(ds, icontrols) * uI(icontrols))
+ bbI(ds)
=e=
x("1",ds);
dely(mesh,ds) ..
x(mesh,ds)
+ ( sum(controls, B(ds, controls) * u(mesh,controls))
+ sum(states, A(ds,states) * x(mesh,states))
+ bb(ds)
) * timespan / N
=e=
x(mesh+1,ds)
+ (xR(ds)) $(ord(mesh) eq N);
model elqp / deluI.uI, delu.u, delx.x, delxR.xR, delv.v,
delvT.vT, delyL.yL, dely.y /;
option limrow = 0;
option limcol = 0;
option solprint = off;
option iterlim = 20000;
option reslim = 14400; /* 4 hours */
file cp /out.pth/;
file qpf /out.mns/;
file cp_obj /obj.pth/;
file qpf_obj /obj.mns/;
scalars
penalty,
object;
uI.l(icontrols) = 0;
u.l(mesh,controls) = 0;
v.l(mesh,duals) = 0;
vT.l(tduals) = 0;
x.l(mesh, states) = 0;
xR.l(states) = 0;
v.l(mesh, duals) = 0;
vT.l(tduals) = 0;
solve elqp using mcp;
penalty = sum (mesh,
sum(duals, v.l(mesh,duals) * (qq(duals)
- sum(controls, D(duals,controls) * u.l(mesh,controls))
- sum(states, C(duals,states) * x.l(mesh,states))
- 0.5 * sum(duals2, Q(duals,duals2) * v.l(mesh,duals2))
) ) ) / N
+ sum(tduals, vT.l(tduals) * (qqT(tduals)
- sum(states, CT(tduals,states) * xR.l(states))
- 0.5 * sum(tduals2, QT(tduals,tduals2) * vT.l(tduals2))
) );
object = sum(mesh,
sum(controls,u.l(mesh,controls) * (pp(controls) +
0.5*sum(controls2,P(controls,controls2)*u.l(mesh,controls2)) ) )
- sum(states, x.l(mesh,states) * cc(states))
) / N
+ sum(icontrols,uI.l(icontrols) * (ppI(icontrols) +
0.5*sum(icontrols2,PI(icontrols,icontrols2)*uI.l(icontrols2))) )
- sum(states, xR.l(states) * ccT(states))
+ penalty;
put cp_obj;
put "ELQP primals:" pcount:5:0 " duals:" dcount:5:0 /;
put "objective function value: " object:15:8 /;
put cp;
$include 'put.pth'
/* start stuff for QP */
putclose cp;
putclose cp_obj;
/* QP formulation described in
R.T. Rockafellar and R.J.-B. Wets,
"A Lagrangian Finite Generation Technique for Solving
Linear-Quadratic Problems in Stochastic Programming",
Math. Prog. Study 28 (1986) 63-93
*/
positive variables
vTdup(tduals),
vTdlo(tduals),
vdup(mesh,duals),
vdlo(mesh,duals);
variables
z, /* objective value */
vv(mesh,duals), /* u'' of RTR paper */
vvT(tduals); /* u'' of RTR paper */
equations
obj,
rfv(mesh,duals),
rfvT(tduals),
rfyL(ds),
rfy(mesh,ds);
obj .. z =e=
sum(icontrols, uI(icontrols) * (ppI(icontrols)
+ 0.5 * sum(icontrols2, PI(icontrols, icontrols2) * uI(icontrols2)) )
)
+ sum(mesh,
sum(controls, u(mesh,controls) * (pp(controls)
+ 0.5 * sum(controls2, P(controls,controls2) * u(mesh,controls2)) ) )
- sum(states, x(mesh,states) * cc(states) )
) / N
- sum(states, xR(states)*ccT(states))
+ sum(mesh,
sum(duals, v_up(duals) * vdup(mesh,duals)
-v_lo(duals) * vdlo(mesh,duals))
)
+ sum(tduals, vT_up(tduals) * vTdup(tduals)
-vT_lo(tduals) * vTdlo(tduals))
+ 0.5 * sum(mesh,
sum(duals, vv(mesh,duals) *
sum(duals2, Q(duals,duals2)*vv(mesh,duals2)))
) / N
+ 0.5 * sum(tduals, vvT(tduals) *
sum(tduals2, QT(tduals,tduals2)*vvT(tduals2)));
rfv(mesh,duals) ..
( sum(controls, D(duals,controls) * u(mesh,controls))
+ sum(states, C(duals,states) * x(mesh,states))
+ sum(duals2, Q(duals, duals2) * vv(mesh,duals2))
- qq(duals)
) / N
+ vdup(mesh,duals) - vdlo(mesh,duals)
=e= 0;
rfvT(tduals) ..
sum(states, CT(tduals,states) * xR(states))
+ sum(tduals2, QT(tduals,tduals2) * vvT(tduals2))
- qqT(tduals)
+ vTdup(tduals) - vTdlo(tduals)
=e= 0;
rfyL(ds) ..
sum(icontrols, BI(ds, icontrols) * uI(icontrols))
+ bbI(ds)
=e=
x("1",ds);
rfy(mesh,ds) ..
x(mesh,ds)
+ ( sum(controls, B(ds, controls) * u(mesh,controls))
+ sum(states, A(ds,states) * x(mesh,states))
+ bb(ds)
) * timespan / N
=e=
x(mesh+1,ds)
+ (xR(ds)) $(ord(mesh) eq N);
model qp / obj, rfv, rfvT, rfyL, rfy /;
uI.l(icontrols) = 0;
u.l(mesh,controls) = 0;
v.l(mesh,duals) = 0;
vT.l(tduals) = 0;
x.l(mesh, states) = 0;
xR.l(states) = 0;
v.l(mesh, duals) = 0;
vT.l(tduals) = 0;
vv.l(mesh,duals) = 0;
vvT.l(tduals) = 0;
vTdup.l(tduals) = 0;
vTdlo.l(tduals) = 0;
vdup.l(mesh,duals) = 0;
vdlo.l(mesh,duals) = 0;
z.l = 999999999;
option iterlim = 60000;
option reslim = 36000; /* 10 hours */
* solve qp using nlp minimizing z;
put qpf_obj;
put "ELQP primals:" pcount:5:0 " duals:" dcount:5:0 /;
put "objective function value: " z.l:15:8 /;
put qpf;
$include 'put.mns'