write dynamic model in python (coming from ampl)

[Hanna Schenk]

Dear all,

I have been using ampl to solve nonlinear optimisation problems with many variables from dynamic models.

in ampl it is easy to write something like this

var x1 {t in 0..T} >= 0

var x2 {t in 0..T} >= 0

var g1 {t in 0..T} = r * x1[t] * (1 - x1[t]/k)

var g2 {t in 0..T} = a * x1[t] * x2[t]

maximise: sum {t in 0..T} ( g2[t] - y[t] )

subject to c1 {t in 0..T}: x2[t+1] = (1 - b) * x2[t] + y[t];
subject to c2 {t in 0..T-1}: x1[t+1] = x1[t] + g1[t] - g2[t];

(a, r, b, k are constants) and solve it with knitro. I have been looking at the python examples and the knitro documentation.

The number of variables in my model depends on T, so I can not write the problem as in the examples/documentation.

Maybe someone can help and tell me whether it is possible (and if so how) to implement this type of problem in python and solve it using knitro

Thanks!

Hanna

[Julien Côté-Massicotte]

Hi Hanna,
Here is an example on how it could be written in Python. It is not complete, but that should give you a good start.
You can pass T, a, r, b and k as userParams to define obj and c.

#*------------------------------------------------------------------*
#*     FUNCTION EVALUATION CALLBACKS                                *
#*------------------------------------------------------------------*


def callbackEvalObj (kc, cb, evalRequest, evalResult, userParams):
    Tvar = userParams.T
 
    if evalRequest.type != KN_RC_EVALFC:
        print ("*** callbackEvalObj incorrectly called with eval type %d" % evalRequest.type)
        return -1
    x = evalRequest.x
   
    for ind in range(Tvar):
        g1[ind] = userParams.r * x[ind] * (1-x[ind]/userParams.k)
        g2[ind] = userParams.a * x[ind] * x[Tvar+ind]
    # Evaluate nonlinear term in objective

    for ind in range(Tvar)
         evalResult.obj += g2[ind] - y[ind]     #multiple ways to write this function
    for ind in range(Tvar-1):
        evalResult.c[ind]= (1-userParams.b)*x[Tvar+ind]+ y[ind] - x[Tvar+1+ind]
        evalResult.c[Tvar+ind]=x[ind] + g1[ind] - g2[ind] - x[ind+1]
    evalResult.c[2*Tvar-1]=... (last constraint for c1, but I'm not sure it exist because x2[t+1] doesn't exist for ind=T)

    return 0

class UserEvalType:
    def __init__ (self, T, a, r, b, k):
        self.T = T
        self.a = a
        self.r = r
        self.b = b
        self.k = k

# Create a new Knitro solver instance.
try:
    kc = KN_new ()
except:
    print ("Failed to find a valid license.")
    quit ()
T=4

a=...

r=...

b=...

k=...
xIndices = KN_add_vars (kc, 2*T)
cIndices = KN_add_cons (kc, 2*T-1)
for x in xIndices:
    KN_set_var_lobnds (kc, x, 0.0)
for c in cIndices:
    KN_set_con_eqbnds (kc, c, 0.0)    #added the left terms to evalResult.c and made it equal to 0.



KN_set_obj_goal (kc, KN_OBJGOAL_MAXIMIZE)

userEval = UserEvalType (T, a, r, b, k)

cbObj = KN_add_eval_callback (kc, evalObj = True, funcCallback = callbackEvalObj)

KN_set_cb_user_params (kc, cbObj, userEval)


Hope this helps!

Best regards,

Julien Côté-Massicotte