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)

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)

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)


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)

T0 <- 3 + cumsum(P0 - P)

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)

n <- 72
d0 <- 3
d1 <- 5
d2 <- 15
d3 <- 100

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)))

require(alabama)
# ?auglag

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
    )
}

solA1 <- auglag(par = P0, fn = funA, hin = hin, control.outer = list(trace = FALSE))
xA1 <- solA1$par
solA1$value

P1 <- solA1$par; T1 <- 3 + cumsum(P1 - P)
plot_gas_solution(T, P, P1, T1)

require(alabama)
solB1 <- auglag(par = P0, fn = funB, hin = hin, control.outer = list(trace = FALSE))
xB1 <- solB1$par
solB1$value

P2 <- solB1$par; T2 <- 5 + cumsum(P2 - P)
plot_gas_solution(T, P, P2, T2)

require(alabama)
solC1 <- auglag(par = P0, fn = funC, hin = hin, control.outer = list(trace = FALSE))
xC1 <- solC1$par
solC1$value

P3 <- solC1$par; T3 <- 3 + cumsum(P3 - P)
plot_gas_solution(T, P, P3, T3)

require(alabama)
solD1 <- auglag(par = P0, fn = funD, hin = hin, control.outer = list(trace = FALSE))
xD1 <- solD1$par
solD1$value

P4 <- solD1$par; T4 <- 3 + cumsum(P4 - P)
plot_gas_solution(T, P, P4, T4)

require(nloptr)
# ?nloptr::auglag

( sol <- nloptr::auglag(x0 = P0, fn = funD, hin = hin, localsolver = "LBFGS") )

P5 <- sol$par; T5 <- 3 + cumsum(P5 - P)
plot_gas_solution(T, P, P5, T5)

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
}

require(DEoptim)
# ?DEoptim

sol <- DEoptim(funDconstr, lower = rep(0, 72), upper = rep(100, 72),
                control=DEoptim.control(trace = FALSE))
sol$optim