Specifying Rate of Change in the Objective Function

Question

I would like to solve some optimization problems with the following objective function

Basically, I want to control the rate of change of x. I am not quite sure what is the most decent way to set up the objective function like this in GEKKO using m.Minimize().

The idea is to approximate the derivative dx/dt using finite difference, let say the most intuitive one: forward difference.

I have tried many different ways, but all of them failed. Below I show all approaches I have done. Note that all the ODEs for x(t) and parameters are specified correctly, and so I do not show all of them here.

(1) The naivest one is something like

obj = 0
for i in range(nt-1):
    obj = obj + ((1/dt)*(x[i+1]-x[i])-0.5)**2

m.Minimize(obj)

The error is: TypeError: 'int' object is not subscriptable.

(2) Then I tried to create a finite difference matrix like (*nt is the number of time points)

q = np.diag(-1*np.ones(nt)) + np.diag(1*np.ones(nt-1),1)
q[-1,-1] = 0
diff = m.Param(value = q)

And then, set the objective to be

obj = ((1/dt)*(diff@x)- 0.5)**2

Basically, I tried to use a similar approach when we optimize the final value, where we define a vector where all entries are zero except the final one. However, I got an error: TypeError: unsupported operand type(s) for @: 'GKParameter' and 'GKVariable'

(3) I read some other previous posts on optimization and so tried something like

q = np.diag(-1*np.ones(nt)) + np.diag(1*np.ones(nt-1),1)
q[-1,-1] = 0

for i in range(nt):
    diff = m.Param(value=q[i])
    m.Minimize(m.abs3((1/dt)*(diff*x)-0.5)**2)

There is no error, but the solver did not even run.

To be honest, I feel like I do not clearly understand how things work inside m.Minimize, even after going through some documentation and previous posts, and so I cannot set the objective function properly as I want. Really appreciate if someone could help elucidate this!

Answer 1

Define a new variable to use a derivative value in the objective function:

y = m.Var(0.5)
m.Equation(y==x.dt())
m.Minimize((y-0.5)**2)

Gekko includes every time point in the horizon to calculate the objective. Below is a simple script that demonstrates how to define the objective function that includes the derivative of x.

from gekko import GEKKO
import numpy as np
m = GEKKO()
m.time = np.linspace(0,10,11)
x = m.Var(0)
y = m.Var(0.5)
m.Equation(y==x.dt())
m.Minimize((y-0.5)**2)

m.options.IMODE=6
m.solve()

import matplotlib.pyplot as plt
plt.figure(figsize=(7,4))
plt.subplot(2,1,1)
plt.plot([0,10],[0.5,0.5],'r-',label='dx/dt=0.5')
plt.plot(m.time,y,'b:',label='y')
plt.legend(); plt.grid()
plt.subplot(2,1,2)
plt.plot([0,10],[0,5],'r-',label='dx/dt=0.5')
plt.plot(m.time,x,'b:',label='x')
plt.legend(); plt.grid()
plt.tight_layout()
plt.show()

You can also adjust the number of NODES in the derivative calculation for higher accuracy with Orthogonal Collocation on Finite Elements. This is set with m.options.NODES=3 between 2 and 6. The default is 2.

It is also possible to switch to IMODE=3 and have more control over the discretization of the differential equation. Here is the equivalent problem with backwards difference Euler's method:

from gekko import GEKKO
import numpy as np
m = GEKKO(remote=False)
n = 11
tm = np.linspace(0,n-1,n)
x = m.Array(m.Var,n)
y = m.Array(m.Var,n)
m.Equation(x[0]==0) # initial conditions
m.Equation(y[0]==0.5)
m.Equations([(tm[i]-tm[i-1])*y[i]==x[i]-x[i-1] 
             for i in range(1,n)])
for i in range(n):
    m.Minimize((y[i]-0.5)**2)

m.options.IMODE=3
m.solve(disp=False)

print(f'y: {y}')
print(f'x: {x}')

Some slight modification is needed for the forward difference formula for approximation of the derivative values:

from gekko import GEKKO
import numpy as np
m = GEKKO(remote=False)
n = 11
tm = np.linspace(0,n-1,n)
x = m.Array(m.Var,n)
y = m.Array(m.Var,n)
m.Equation(x[0]==0) # initial condition
m.Equations([(tm[i+1]-tm[i])*y[i]==x[i+1]-x[i] 
             for i in range(0,n-1)])
m.Equation(y[-1]==y[-2]) # final condition
for i in range(n):
    m.Minimize((y[i]-0.5)**2)

m.options.IMODE=3
m.solve(disp=False)

print(f'y: {y}')
print(f'x: {x}')

They both produce the same answer as above:

y: [[0.5] [0.5] [0.5] [0.5] [0.5] [0.5] [0.5] [0.5] [0.5] [0.5] [0.5]]
x: [[0.0] [0.5] [1.0] [1.5] [2.0] [2.5] [3.0] [3.5] [4.0] [4.5] [5.0]]

where the derivative of x is y=0.5 and x is a line that starts at zero and has a slope of 0.5.