An Example in Constraint Optimization

Introduction

The task is to manage energy consumption for a local site, for example a neighborhood of households or an industrial area. The predicted gas consumption for the next 72 hours (3 days) is given in the following vector, say in cubic meters or similar.

In [1]:
T <- 1:72
P <- c( 96.6, 96, 98.4, 99.1, 98, 100.6, 100.7, 99.4, 99.5, 98.9, 98, 98.3, 96.3, 
        98.4, 95.5, 95.3, 95.8, 94.6, 94.4, 93.8, 95.2, 94, 96.6, 98.2, 96.7, 99.2, 
        100.6, 98.6, 99.4, 102, 101.3, 102.5, 102.5, 102.1, 102.2, 100.9, 100.8, 
        97.6, 98.1, 97.4, 95.4, 96.9, 95.3, 97.7, 95.5, 97, 96.9, 96, 95.3, 98, 
        97.6, 99.7, 99.9, 101.2, 100, 98.6, 99.2, 97.4, 97.7, 97.1, 96.5, 95.5, 
        96.4, 93, 95, 95.4, 95.3, 95.4, 95.3, 96.4, 96.9, 96.7)

The following plot displays this prediction as a step curve. The local supplier of energy buys gas from a global supplier. The contract says he will get up to 100 [m^3] per hour to a fixed price, if more gas is consumed the price will go up significantly. The red line marks this upper limit.

In [2]:
plot(c(0, 72), c(90, 105), type = "n",
    main = "Predicted gas consumption", xlab = "Time", ylab = "Gas")
grid(); abline(h = 100, col = "red")
lines(T, P, col = "green", type = "s", lwd = 2)

To avoid these extra costs the local provider will employ some kind of gas storage, e.g. a gas tank (or the so-called "net aspiration"). When gas consumption in the area is low, the tank can be filled, and when consumption is greater than 100 the gas storage will be used to contribute the missing amount.

Of course, the gas storage facility has some technical constraints that need to be taken into account. The maximum amount of gas to be stored in the tank shall be 10 [m^3], so the actual amount of gas in the tank will vary between 0 and 10. Also, the amount of gas that will be pumped in or taken out has be less or equal to 5 [m^3] per hour.

Imagine the planned amount of gas bought from the global supplier is

In [3]:
P0 <- c(
100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0,  99.2,  98.0,  98.3,
 96.3,  98.4,  95.5,  95.3,  95.8,  94.6,  94.4,  93.8,  95.2,  94.0,  96.6,  98.2,
 96.7,  99.2, 100.0,  99.2,  99.4, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0,
100.0, 100.0, 100.0, 100.0, 100.0,  99.7,  95.3,  97.7,  95.5,  97.0,  96.9,  96.0,
 95.3,  98.0,  97.6,  99.7,  99.9, 100.0, 100.0,  99.8,  99.2,  97.4,  97.7,  97.1,
 96.5,  95.5,  96.4,  93.0,  95.0,  95.4,  95.3,  95.4,  95.3,  96.4,  96.9,  96.7)

then the combined display of predicted and planned gas consumption is shown in the following plot:

In [4]:
plot(c(0, 72), c(90, 105), type = "n",
    main = "Predicted gas consumption", xlab = "Time", ylab = "Gas")
grid(); abline(h = 100, col = "red")
lines(T, P, col = "green", type = "s", lwd = 2)
lines(T, P0, col = "navy", type = "s", lwd = 2, lty = 2)

Can this plan be satisfied by utilizing the gas storage available to the local supplier? Assuming the gas tank is half-filled with 5 [m^3], what will be the level in the tank during the planned 72 hours? Obviously, the filling level is calculated as 5 + cumsum(P0 - P). So with this planning the tank level varies as follows:

In [5]:
T0 <- 3 + cumsum(P0 - P)

From now on we will plot the (predicted or planned) gas consumption and the level in the gas tank -- i.e. the volume of stored gas -- side by side in one graph. For this greedy approach the figure looks as follows:

In [6]:
plot_gas_solution <- function(T, P, P0, T0) {
    par(mfrow = c(2, 1), ann = FALSE, mai = c(0.5, 0.5, 0.25, 0.5))
    plot(T, P, type = "s", col = "green", ylim = c(90, 105)); grid()
    lines(T, P0, type = "s", col = "blue")
    plot(T, T0, type = "s", col = "red", ylim = c(0, 20)); grid()
}

plot_gas_solution(T, P, P0, T0)

We see that the gas tank is completely filled most of the time and that it can handle the problem of a missing gas volume during the 30--40 hours ahead.

Problem Formulation

The task now is to formulate this as an optimization problem such that, when given a prediction of gas consumption for 72 hours, a purchase plan for buying gas from the global provider will be generated. The plan shall make it possible to never buy more than a certain amount of gas, and the gas tank shall be used in a moderate way.

Decision Variables

Let $n$ be the number of time periods, $p_i, i=1, ...,72$ be the predicted consumption values, and $q_i$ the unknown planning values to be optimized during the optimization procedure.

How can a feasible starting point x0 be found? The same way P0 above was found. That is, immediately fill up the tank and keep it filled as much as possible. If this leads to a non-feasible solution -- because the storage tank is not big enough -- than we know that the whole optimization problem is not solvable. It is a very nice feature of this optimization problem that feasibility can so easily be determined, and if applicable a feasible starting point is constructed.

Constraints

The following constants need to be defined:

  • $n$, the number of predicted resp. planned time periods.

  • $d_0$: initial level of the gas tank.

  • $d_1$: maximum level change per time unit.

  • $d_2$: maximum level of gas tank (minimum is 0).

  • $d_3$: maximum allowed amount of gas from the global supplier.

In [7]:
n <- 72
d0 <- 3
d1 <- 5
d2 <- 15
d3 <- 100

With these values the constraints are defined as

  1. maximum supply to be ordered from the global supplier:
    $$ q_i \leq d_3 \qquad \mathrm{for\ all} \; i = 1 \ldots n $$

  2. restricted amount to be taken from or filled into gas tank: $$ |q_i - p_i| \le d_1 \qquad \mathrm{for\ all} \; i = 1 \ldots n$$

  3. level in storage tank restricted by capacity: $$ 0 \le d_0 + \sum_{k=1}^{j} (q_k - p_k) \leq d_2 \qquad \mathrm{for\ all} \; j = 1 \ldots n$$

Objective

The objective function will describe the overall goal that is to be reached.
Four different goals will be considered and compared:

  • (A) Minimize the sum of differences between predicted and planned purchase of gas from the global supplier: $$ \min ! \sum_{j=1}^n |p_j - q_j| $$

  • (B) Minimize the sum of squares of these differences: $$ \min ! \sum_{j=1}^n (p_j - q_j)^2 $$ (Larger differences will have a larger influence on the objective.)

  • (C) Minimize the difference between subsequent supplies: $$ \min ! \sum_{j=1}^{n-1} |q_{j+1} - q_{j}| $$

  • (D) Minimize usage of the gas tank as much as possible: $$ \min ! \sum_{k=1}^{n} | \sum_{j=1}^{k} (p_j - q_j)|$$

We will see that very different solutions will be returned with each of these objectives. Of course, it is up to the customer of such an optimization package to decide which solution and thus which objective function he prefers.

The objective functions will be the same for all solvers, so we will define them here.

In [8]:
funA <- function(q) sum(abs(P - q))

funB <- function(q) sum((P - q)^2)

funC <- function(q) sum(abs(diff(q)))

funD <- function(q) sum(abs(cumsum(P - q)))

Nonlinear Optimization Solvers in R

The function constrOptim() in base R is a solver for nonlinear optimization problems with linear constraints. It implements an adaptive logarithmic barrier algorithm, using optim() as the underlying solver. Its usage is as follows:

constrOptim(theta, f, grad, ui, ci, mu = 1e-04, control = list(),  
                method = if (is.null(grad)) "Nelder-Mead" else "BFGS", 
                outer.iterations = 100, outer.eps = 1e-05, ..., hessian = FALSE)

But our constraints are nonlinear, so we need a more powerful tool. The appropriate technique here is the augmented Lagrangian. Solvers available in R packages that apply this approach are

  • nloptr
  • alabama
  • Rsolnp
  • Rdonlp2

Unfortunately, different optimization solvers will require to receive the the constraints and objective function in different ways. The following solvers for non-linear constraint optimization will be tried on the problem: 'alabama', 'nloptr'.

Because all these solvers are local solvers, i.e. will stop in local minima, we will also try a global (or stochastic) solver and compare it with solutions found with local solvers.

The alabama package

Solver auglag() in package 'alabama' optimizes smooth nonlinear objective functions with nonlinear constraints. Equality and inequality constraints are allowed. As the problem formulated above does not involve equality constraints, the full power of an augmented Lagrangian approach is not required.

The general call to auglag() looks like this:

auglag(par, fn, gr, hin, hin.jac, heq, heq.jac, 
       control.outer=list(), control.optim = list(), ...)

with the following parameters

Argument Meaning
par initial vector of variable values
fn nonlinear objective function
gr gradient of the objective function
hin function specifying the inequality constraints
... ...

And control.optim a list of control parameters, the same as those used in optim(). The default method is "BFGS". At the moment we rely on numerical gradients and jacobians. For more information see ?auglag.

In [9]:
require(alabama)
# ?auglag
Loading required package: alabama
Loading required package: numDeriv

The constraint inequalities need to be defined through hin[j] >= 0 for all j. Therefore, we define function hin() as

In [10]:
hin <- function(q) {
    c(d3 - q,                   # q[i] <= d3
      d1 - abs(P - q),          # abs(P - q) <= d1
      d0 + cumsum(q - P),       # 0 <= d0 + cumsum(q - P)
      d2 - d0 - cumsum(q - P)   # d0 + cumsum(q - P) <= d2
    )
}

Objective funA

In [11]:
solA1 <- auglag(par = P0, fn = funA, hin = hin, control.outer = list(trace = FALSE))
xA1 <- solA1$par
solA1$value
Out[11]:
31.8000008650291

Display the solution and the level in the storage tank.

In [12]:
P1 <- solA1$par; T1 <- 3 + cumsum(P1 - P)
plot_gas_solution(T, P, P1, T1)

In this solution, the level in the gas tank is kept at a minimum and only raised when it is absolutely necessary to provide for peaks in the predicted consumtion. Compare this with a solution to problem (B).

Objective funB

In [13]:
require(alabama)
solB1 <- auglag(par = P0, fn = funB, hin = hin, control.outer = list(trace = FALSE))
xB1 <- solB1$par
solB1$value
Out[13]:
38.3316001342559
In [14]:
P2 <- solB1$par; T2 <- 5 + cumsum(P2 - P)
plot_gas_solution(T, P, P2, T2)

In this solution, the level in the gas tank is only raised slowly and only when it is absolutely necessary to provide for peaks in the predicted consumption.

Objective funC

For problem type (C) it is the goal to minimize the differences between gas intakes from the global supplier.

In [15]:
require(alabama)
solC1 <- auglag(par = P0, fn = funC, hin = hin, control.outer = list(trace = FALSE))
xC1 <- solC1$par
solC1$value
Out[15]:
8.13892864328147
In [16]:
P3 <- solC1$par; T3 <- 3 + cumsum(P3 - P)
plot_gas_solution(T, P, P3, T3)

In this case the level in the gas tank varies quite much.

Objective funD

For problem type (D) it is the goal to minimize the usage of the gas tank as much as possible.

In [17]:
require(alabama)
solD1 <- auglag(par = P0, fn = funD, hin = hin, control.outer = list(trace = FALSE))
xD1 <- solD1$par
solD1$value
Out[17]:
98.2018402571746
In [18]:
P4 <- solD1$par; T4 <- 3 + cumsum(P4 - P)
plot_gas_solution(T, P, P4, T4)

Interestingly, in this case the level in the gas tank is kept to the level it had at the start, is raised if there is a need ahead, and then kept constant at the initial value.

The nloptr package

In [19]:
require(nloptr)
# ?nloptr::auglag
Loading required package: nloptr

Attaching package: ‘nloptr’

The following object is masked from ‘package:alabama’:

    auglag

In [20]:
( sol <- nloptr::auglag(x0 = P0, fn = funD, hin = hin, localsolver = "LBFGS") )
Out[20]:
$par
  1. 96.5989061873907
  2. 96.0011077076734
  3. 98.4005555051972
  4. 99.0994435528671
  5. 98.6457084631035
  6. 99.9999665185044
  7. 99.9989653881265
  8. 99.9999945813929
  9. 99.5583919020994
  10. 98.8972276331535
  11. 98.0009020987691
  12. 98.2961273810968
  13. 96.3028340331549
  14. 98.3994499819561
  15. 95.4986924879238
  16. 95.3021315721442
  17. 95.8031115236202
  18. 94.5945708108803
  19. 94.3826359706308
  20. 93.8171263893068
  21. 95.5554938941084
  22. 95.177057627338
  23. 99.0708686195802
  24. 99.9993825930707
  25. 99.9995861710969
  26. 99.9998976741525
  27. 99.9999035669954
  28. 99.9999965355193
  29. 99.9999971797784
  30. 99.9999968280422
  31. 99.9999968335165
  32. 99.9999971900111
  33. 99.9999975794571
  34. 99.9999983198966
  35. 99.9999975109418
  36. 99.9999972771428
  37. 99.9999970953074
  38. 99.9995262961189
  39. 98.6998496346533
  40. 97.4003692905741
  41. 95.4009031078898
  42. 97.0947559853778
  43. 95.1032249110913
  44. 97.7023092405544
  45. 95.4986695177783
  46. 96.9985083373798
  47. 96.9038493035109
  48. 95.9981196938435
  49. 95.3013725692337
  50. 97.9983533628279
  51. 97.6000359091605
  52. 99.9906361085301
  53. 99.9999266231285
  54. 99.9999269022413
  55. 99.9999766305059
  56. 99.4096328081496
  57. 99.1992908615297
  58. 97.4009138506745
  59. 97.7004012144889
  60. 97.0994298067916
  61. 96.5001398156122
  62. 95.5024442058138
  63. 96.3975353784967
  64. 92.9988157198024
  65. 95.0014711921282
  66. 95.3986000873134
  67. 95.2997130317036
  68. 95.401664824706
  69. 95.3009026931983
  70. 96.3985619407515
  71. 96.9010691876763
  72. 96.6990774390717
$value
99.743653222224
$iter
729
$global_solver
'NLOPT_LD_AUGLAG'
$local_solver
'NLOPT_LD_LBFGS'
$convergence
4
$message
'NLOPT_XTOL_REACHED: Optimization stopped because xtol_rel or xtol_abs (above) was reached.'
In [21]:
P5 <- sol$par; T5 <- 3 + cumsum(P5 - P)
plot_gas_solution(T, P, P5, T5)

This agrees quite exactly with the solution found by alabama::auglag !

Some Remarks

Global Solvers

There are several global (or: stochastic) solvers in R packages that try to avoid falling into or getting stuck in local minima. There is absolutely no guaranty they will find the global minimum. The problem here with 72 dimensions may anyway be to demanding.

For realizing constraints for application in global solvers the objective function has to be changed to include the constraints.

In [22]:
funDconstr <- function(q) {
    if (# any(q > d3) ||
        any(abs(P - q) > d1 + 0.001) ||
        any(d0 + cumsum(q - P) < -0.001) || 
        any(d0 + cumsum(q - P) > d2 + 0.001) ){
        v <- Inf
    } else {
        v <- sum(abs(cumsum(P - q)))
    }
    v
}
In [23]:
require(DEoptim)
# ?DEoptim
Loading required package: DEoptim

DEoptim package
Differential Evolution algorithm in R
Authors: D. Ardia, K. Mullen, B. Peterson and J. Ulrich

In [24]:
sol <- DEoptim(funDconstr, lower = rep(0, 72), upper = rep(100, 72),
                control=DEoptim.control(trace = FALSE))
sol$optim
Out[24]:
$bestmem
par1
13.0816068990255
par2
37.7001394098625
par3
58.0040045467836
par4
16.1892402737211
par5
59.4149708747864
par6
38.4608807369184
par7
46.5962492162362
par8
78.7978068853282
par9
38.4119799384688
par10
72.4666442910566
par11
92.4671353143437
par12
37.0134336412961
par13
80.599758140852
par14
84.7988849632056
par15
36.2438437682034
par16
75.8663427920953
par17
47.8722647763789
par18
37.7616989056118
par19
40.7056774892468
par20
60.5002468574715
par21
14.2377811267308
par22
22.0176759408787
par23
30.7340047310155
par24
36.8046016898006
par25
37.8502323596373
par26
85.9464909826515
par27
61.7429656586262
par28
81.8412055269053
par29
57.7491433126852
par30
49.7165983660809
par31
20.1104662033919
par32
29.1323970043018
par33
72.1030247409978
par34
13.3235943363391
par35
58.1106220516493
par36
54.9491913483415
par37
66.3247101630556
par38
80.9321391788396
par39
61.9171395651384
par40
53.3641762266022
par41
70.7569820724694
par42
18.6819009017199
par43
39.8277879252831
par44
74.5165301711985
par45
48.2873672686928
par46
50.7920337592744
par47
49.758368138106
par48
78.4112193852888
par49
72.5511527016619
par50
4.26100119837755
par51
16.7142282147089
par52
19.2892339080572
par53
0.228710095398128
par54
60.4324849259658
par55
92.1863871728705
par56
44.0147380810231
par57
97.0953421540852
par58
82.3845442345574
par59
67.2124875378484
par60
99.3674361171424
par61
85.4644238948822
par62
88.6985497146457
par63
49.2805518675596
par64
62.1778574073687
par65
55.3358704800553
par66
58.4999817046718
par67
45.6047700411342
par68
45.6557461931735
par69
12.9174374365999
par70
66.2025419296697
par71
73.6915599992273
par72
89.7816472982302
$bestval
Inf
$nfeval
402
$iter
200

We can see that global solvers will have severe problems with constraints if the feasible points are rare in the solution space, or if the best solutions lie on or near the boundary.

MATLAB Solution

In MATLAB, the standard solver for constrained nonlinear optimization (i.e., 'nonlinear programming') is fmincon from the Optimization Tolbox. It incorporates active-set, sqp, trust-region, and interior-point algorithms.

The solver fmincon is called with the following syntax and parameters:

fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon)

where obj is the objective function, x0 the starting point, A*x <= b the linear inequality and Aeq*x = beq the linear equality constraints, lb and ub the lower and upper bounds with lb <= x <= ub.
nonlcon defines the nonlinear equality and inequality constraints c(x) <= 0 and ceq(x) = 0 by returning two vectors c and ceq: [c,ceq] = nonlcon(x).

The user interface for a demo system to discuss with engineers and customers was developed utilizing MATLAB's GUI editor and is shown in the following figure.

MATLAB solution

One can see in this figure that fmincon finds the same solutions as auglag in R above.

A point to clarify will be what happens when we start with other feasible points.

Modeling Language

R and MATLAB do not have optimization modeling languages that would enable the user to formulate optimization problems independently of what the API of different solvers will require. Commercial programs such as AMPL, GAMS, or Gurobi each have their own (algebraic) modeling languages. Goals for applying a modeling language are, e.g., user friendliness, mathematical notation, solver independence, and getting access to advanced optimization techniques.

The following is an AMPL model for our minimization problem with target function funB, but extended with a request that the gas tank should be filled up to a level of 5.

AMPL Model

param n > 0 integer;    # n = 72 time units
param p  {i in 1..n};   # predicted consumption

param d0;   # initial gas tank storage
param d1;   # max level change per time unit
param d2;   # max level of gas tank
param d3;   # max amount of gas from supplier

var x {i in 1..n} >= 0; # solution

minimize target:
    sum {j in 1..n} (p[j] - x[j]) * (p[j] - x[j]);

subject to rule1 {i in 1..n}: x[i] <= d3;
subject to rule2 {i in 1..n}: abs(x[i] - p[i]) <= d1;
subject to rule3 {j in 1..n}:
    sum {k in 1..j} (x[k] - p[k]) <= d2 - d0;
subject to rule4 {j in 1..n}:
    sum {k in 1..j} (x[k] - p[k]) >= -d0;

AMPL Data

param n := 72;

param d0 := 3;
param d1 := 5;
param d2 := 15;
param d3 := 100;

param: p :=
    1  96.6  2  96.0 ... 72  96.7;

AMPL Commands

ampl: option solver minos;
ampl: option minos_options 'iterations_limit=5000 superbasics_limit=100';

ampl: model gas_ex.mod;
ampl: data gas_ex.dat;

ampl: solve;
ampl: display x;

Sending this to AMPL will call the MINOS solver on the model and data and will result in the following output:

MINOS 5.51: iterations_limit=5000
superbasics_limit=100
MINOS 5.51: optimal solution found.
237 iterations, objective 38.3316
Nonlin evals: obj = 352, grad = 351, constrs = 352, Jac = 351.

x [*] :=
 1  97.108   13  96.808   25  97.208   37 100       49  95.375   61  96.5
 2  96.508   14  98.908   26  99.708   38  97.675   50  98.075   62  95.5
 3  98.908   15  96.008   27 100       39  98.175   51  97.675   63  96.4
 4  99.608   16  95.808   28  99.108   40  97.475   52  99.775   64  93
 5  98.508   17  96.308   29  99.908   41  95.475   53  99.975   65  95
 6 100       18  95.108   30 100       42  96.975   54 100       66  95.4
 7 100       19  94.908   31 100       43  95.375   55 100       67  95.3
 8  99.908   20  94.308   32 100       44  97.775   56  98.6     68  95.4
 9 100       21  95.708   33 100       45  95.575   57  99.2     69  95.3
10  99.408   22  94.508   34 100       46  97.075   58  97.4     70  96.4
11  98.508   23  97.108   35 100       47  96.975   59  97.7     71  96.9
12  98.808   24  98.708   36 100       48  96.075   60  97.1     72  96.7
;

This corresponds to the solution P2 that we found above.

rneos: XML-RPC Interface to NEOS

"Within this package the XML-RPC API to NEOS is implemented. This enables the user to pass optimization problems to NEOS and retrieve results within R."

The NEOS Server is an internet-based service for solving numerical optimization problems. It provides free access to more than 60 state-of-the-art (commercial and non-commercial) solvers. Optimization problems need to be formulated in AMPL or GAMS syntax. Results will be returned as Web pages.

With the above model and data as files gas_ex.mod and gas_ex.dat, an interaction with the CONOPT solver on NEOS would look like:

require(rneos)
# geting a template for category and solver
temple <-NgetSolverTemplate(category = "nco", solvername = "CONOPT",
                            inputMethod = "AMPL")
# setting model and data file
modf <- "gas_ex.mod"; datf <- "gas_ex.dat"
modc <- paste(paste(readLines(modf), collapse = "\n"), "\n")
datc <- paste(paste(readLines(datf), collapse = "\n"), "\n")

# create list object
argslist <- list(model = modc, data = datc,
                 commands = "", comments = "Gas example")
## create XML string
xmlstring <- CreateXmlString(neosxml = temple, cdatalist = argslist)

## submit job to the NEOS solver
neosjob <- NsubmitJob(xmlstring, user = "hwb", interface = "gas_ex",
                      id = 8237, nc = CreateNeosComm())
neosjob
# The job number is: 3838832 
# The pass word is : wBgHomLT 

# getting info about job
NgetJobInfo(neosjob)            # "nco"   "MINOS" "AMPL"  "Done" 
NgetFinalResults(neosjob)

The available NEOS solvers and the required modeling language(s) can be found on the following NEOS solvers page http://www.neos-server.org/neos/solvers/ at the University of Wisconsin Institutes for Discovery.