Using a variable dependent function in constraints

Viewed 88

I'm trying to optimize the following problem with gekko: enter image description here

The problem is that the dynamics of the problem contain a term, kappa, that is dependent on s but the only relationship I am provided between them are data of kappa at every discrete points in s.

I'm also trying to formulate my problem in both discrete and continuous form to see which one runs faster. My attempt is to estimate kappa as a function of s using cspline but I'm having errors implementing it to the constraints in both. You can see the code below:

In discrete time:

now = time.time()
p = GEKKO(remote=False)

sg = []
rhog = []
thetag = []
betag = []
vg = []
ag = []
curvg = [] #kappa

#parameters
T = p.Const(value=0.5)
vref = p.Const(value=3)
nT = 11

#kappa and s reference
arcref = refpath[2,minid:]-refpath[2,minid]
curvref = refpath[3,minid:]

for i in range(nT):
    #Set variables
    sg.append(p.Var(value = sinit))
    rhog.append(p.Var(value = rhoinit))
    thetag.append(p.Var(value = 0))
    betag.append(p.Var(value = 0))
    vg.append(p.Var(value = 0))
    ag.append(p.Var(value = 0))
    curvg.append(p.Var())
    
    #Set objective function
    p.Obj(rhog[i]**2 + (vg[i]- vref)**2 + ag[i]**2)
    
    #Estimate kappa
    p.cspline(sg[i], curvg[i], arcref, curvref)

for i in range(nT-1):
    #Set constraints
    p.Equation( vg[i+1] == vg[i] + T*ag[i] )
    p.Equation( sg[i+1] == sg[i] + T*vg[i]*p.sin(thetag[i]) )
    p.Equation( rhog[i+1] == rhog[i] + T*vg[i]*p.cos(thetag[i])/(1-rhog[i]*curvg) )
    p.Equation( thetag[i+1]-betag[i+1]+curvg[i+1]*sg[i+1] == thetag[i]-betag[i]+curvg[i]*sg[i] )

p.options.IMODE = 3
p.options.SOLVER = 3
p.options.WEB = 0

p.solve(disp=False, debug=True)
print(f'run time: {time.time()-now}')

and I got the error

Exception: @error: Model Expression
 *** Error in syntax of function string: Missing operator

Position: 2                   
 0,0,0,0,0,0,0,0,0,0,0
  ?

In continuous time:

now = time.time()
p = GEKKO(remote=False)

#Set variables
sg = p.Var(value=sinit)
rhog = p.Var(value=rhoinit)
vg = p.Var(value=0)
ag = p.Var(value=0)
thetag = p.Var(value=0)
betag = p.Var(value=0)

curvg = p.Var()
#Parameters
T = 0.5
vref = p.Const(value=3)
nT = 11

#kappa and s reference
arcref = refpath[2,minid:]-refpath[2,minid]
curvref = refpath[3,minid:]

p.time = np.linspace(0,nT,nT*2)

p.options.IMODE = 6
p.options.SOLVER = 3
p.options.WEB = 0

p.Obj( rhog**2 + (vg - vref)**2 + ag**2 )

#set kappa estimate
p.cspline(sg, curvg, arcref, curvref)

p.Equation(sg.dt()==vg*p.cos(thetag)/(1-rhog*curvg))
p.Equation(rhog.dt()==vg*p.sin(thetag))
p.Equation(vg.dt()==ag)
p.Equation(thetag.dt()==betag.dt()-curvg*sg.dt())

p.solve(disp=False, debug=True)
print(f'run time: {time.time()-now}')

and I got the error:

Error: free(): invalid size

Program received signal SIGABRT: Process abort signal.

Backtrace for this error:
#0  0x6c278f
#1  0x6aacd0
#2  0x7f95ba478fcf
#3  0x7f95ba478f47
#4  0x7f95ba47a8b0
#5  0x7f95ba4c3906
#6  0x7f95ba4ca979
#7  0x7f95ba4d1e9b
#8  0x4f54a6
#9  0x4f66cb
#10  0x50802b
#11  0x5165af
#12  0x52661f
#13  0x449776
#14  0x449d6f
#15  0x4472e2
#16  0x66d119
#17  0x4026ec
#18  0x7f95ba45bb96
#19  0x40275c
#20  0xffffffffffffffff

Error: 'results.json' not found. Check above for additional error details
1 Answers

I couldn't replicate your problem because it was missing some of the inputs. Here is a version that runs successfully in both Windows and on Linux.

from gekko import GEKKO
import time
import numpy as np

now = time.time()
p = GEKKO(remote=False)

#Set variables
sg = p.Var()
rhog = p.Var()
vg = p.Var(value=0)
ag = p.Var(value=0)
thetag = p.Var(value=0)
betag = p.Var(value=0)

curvg = p.Var()
#Parameters
T = 0.5
vref = p.Const(value=3)
nT = 11

#kappa and s reference
arcref = [-1,0,1,2,3,4,5,6]
curvref = [1,2,1,2,1,2,1,2]

p.time = np.linspace(0,nT,nT*2)

p.options.IMODE = 6
p.options.SOLVER = 3
p.options.WEB = 0

p.Obj( rhog**2 + (vg - vref)**2 + ag**2 )

#set kappa estimate
p.cspline(sg, curvg, arcref, curvref)

p.Equation(sg.dt()==vg*p.cos(thetag)/(1-rhog*curvg))
p.Equation(rhog.dt()==vg*p.sin(thetag))
p.Equation(vg.dt()==ag)
p.Equation(thetag.dt()==betag.dt()-curvg*sg.dt())

p.solve(disp=False, debug=True)
print('sg: ' + str(sg.value))
print('curvg: ' + str(curvg.value))
print(f'run time: {time.time()-now}')

I changed the cspline inputs to sample data and initialized sg and rhog to defaults (0). Make sure that the cspline has a monotonically increasing dependent variable.

Related