gekko dynamic optimization - convergence problem

[Ricardo LUNA HERNANDEZ]

Dear Jonh and Logan Beal,

I have had problem with the converge of dynamic optimization in a rice drying model (attached the code).

I have solved this problem in AMPL before using a IPOPT solver and obtained good results. However, with Gekko does not obtained convergence yet. My idea is compare both softwares results.

Then, i would like your help to solve this problem and identify the error.

Finally, i have a question:

Where can see the output specifications of the solver (the numbers 0,-1 or 1)?, for example:

Thanks in advance

Best regards,

Ricardo Luna

[John Hedengren]

Dear Ricardo,

I couldn't get IPOPT to work but APOPT was able to converge. Let me know if you have any questions about the model. I also included a small DCOST on the MVs. This smooths out the MV profiles and doesn't make movements unless there is a reason (decrease the objective function sufficiently). Other changes are in bold. I added an equation for xm because you had only initialized that value but not included an equation. I also changed xm==0.13 to xm<=0.13 and added a soft constraint (1000 * (xm-0.13)**2 * final) to help the optimizer find a solution.

Without DCOST:

With DCOST (1e-5):

# -*- coding: utf-8 -*-

"""

Created on Wed Sep 26 18:18:45 2018

 

@author: Ricardo

"""

 

from gekko import GEKKO

import matplotlib.pyplot as plt

import numpy as np

 

# Initialize Model

m = GEKKO()

#help(m)

 

nt = 21;

tf = 2.0; # hrs

m.time = np.linspace(0,tf,nt)

 

# parameters

 

B10 = m.Param(value = 0.01316*3600);        # (kg dry matter)* m^-3*s^-1

B11 = m.Param(value = 0.3083);              # (kg water)^-1*(kg dry matter)*°C^-1

 

B20 = m.Param(value = 2.304e-09*3600);

B21 = m.Param(value = 0.0442);

 

C1  = m.Param(value = 0.319)

C2  = m.Param(value = 0.0493)

C3  = m.Param(value = 1.8994)

C4  = m.Param(value = 2.5457)

C5  = m.Param(value = 65)

 

Cpg = m.Param(value = 1300)

Cpw = m.Param(value = 4210)

 

Ea = m.Param(value = 1.657e+05)

K0 = m.Param(value = 1.56e+27*3600)

 

Lv = m.Param(value = (2608.8-251.2)*1e+03)

R = m.Param(value = 8.32)

rog = m.Param(value = 1500)

Ssg = m.Param(value = 2000)

Tau1 = m.Param(value = 0.6)

Tau2 = m.Param(value = 0.4)

 

E = m.Param(value = 16.5362)

F = m.Param(value = 3985.44) # antoine parameters

G = m.Param(value = -38.9974)

 

# manipulates variable

Ta = m.MV(value = 60.0, lb = 40.0, ub= 80.0)

Ta.STATUS=1;

Ta.DCOST=1e-5;

Hr = m.MV(value = 0.5,  lb = 0.05, ub = 0.8)

Hr.STATUS=1;

Hr.DCOST=1e-5;

 

# states variables

x1 = m.Var(value = 0.27,lb = 0)

x2 = m.Var(value = 0.27,lb = 0)

Tg = m.Var(value = 20.0,lb = 0)

q  = m.Var(value = 80.0,lb = 0, ub = 100)

xm = m.Var(value = x1*Tau1+x2*Tau2, lb = 0)

 

# constitutive equations

K     = m.Intermediate(-K0*(x2-x1)**5*m.exp(-Ea/(R*(Tg+273.15))))

 

beta1 = m.Intermediate(B10*m.exp(B11*xm*Tg))

beta2 = m.Intermediate(B20*m.exp(B21*Ta))

alpha = m.Intermediate(C5*Lv*beta2)

Aw    = m.Intermediate(m.exp(-m.exp((C1-x2)/C2)/(C3*(Tg-C4))))

Pgsat = m.Intermediate(m.exp(E-F/(G+(Tg+273.15)))*1e+03)

pg    = m.Intermediate(Pgsat*Aw)

pa    = m.Intermediate(Hr*m.exp(E-F/(G+(Ta+273.15)))*1e+03)

 

# differential equations

m.Equation(x1.dt()*(rog*Tau1) == beta1*(x2-x1))

m.Equation(x2.dt()*(rog*Tau2**2) == (beta2*Ssg)*(pa-pg)-x1.dt()*Tau1)

m.Equation(Tg.dt()*(rog*(Cpg+Cpw*xm)) == (alpha*Ssg*(Ta-Tg)+beta2*Ssg*(pa-pg)*Lv))

m.Equation(q.dt()  == -K*q**2)

m.Equation(xm==x1*Tau1+x2*Tau2)

 

# final point

p = np.zeros(nt)

p[-1]=1.0;

final = m.Param(value=p)

 

# equality constraint

m.Equation(xm*final <= 0.13)

# inequality constraint

m.Equation(q.dt() <= 0)

 

# objective function

m.Obj(-final*q)

m.Obj(1000 * (xm-0.13)**2 * final)

m.options.IMODE = 6;  # dynamic optimization

m.options.NODES = 3;  # number of colocations points

m.options.SOLVER = 1; # 1: APOPT, 2:BPOPT, 3:IPOPT

m.options.MAX_ITER = 300 # change maximum iterations

 

m.solve()

 

 

#print(xm.value)

 

plt.figure(figsize=(12,5))

plt.subplot(3,2,1)

plt.plot(m.time,Ta.value,'k:',LineWidth=2,label=r'$Ta$')

plt.xlabel('Time (h)')

plt.ylabel('$Ta$ (°C)')

 

plt.subplot(3,2,2)

plt.plot(m.time,Hr.value,'r:',LineWidth=2,label=r'$Hr$')

plt.xlabel('Time (h)')

plt.ylabel('$Hr$ (°C)')

 

 

plt.subplot(3,2,3)

plt.plot(m.time,Tg.value,'k:',LineWidth=2,label=r'$Tg$')

plt.xlabel('Time (h)')

plt.ylabel('$Tg$ (°C)')

 

plt.subplot(3,2,4)

plt.plot(m.time,q.value,'r:',LineWidth=2,label=r'$q$')

plt.xlabel('Time (h)')

plt.ylabel('$q$ (%)')

 

plt.subplot(3,2,5)

plt.plot(m.time,x1.value,'k:',LineWidth=2,label=r'$x1$')

plt.plot(m.time,x2.value,'r:',LineWidth=2,label=r'$x2$')

plt.legend(loc='best')

plt.xlabel('Time (h)')

plt.ylabel('$x$')

 

plt.subplot(3,2,6)

plt.plot(m.time,xm.value,'g-',LineWidth=2,label=r'$x_m$')

plt.plot(m.time[-1],[0.13],'bo',label='end constraint')

plt.legend(loc='best')

plt.xlabel('Time (h)')

plt.ylabel('$xm$')

 

plt.show()

 

You can compare the difference between the two cases and the DCOST doesn't make the solution that much worse. Nice application!

Best regards,

John Hedengren

--
--
APMonitor user's group e-mail list.
- To post a message, send email to apmo...@googlegroups.com
- To unsubscribe, send email to apmonitor+...@googlegroups.com
- Visit this group at http://groups.google.com/group/apmonitor

---
You received this message because you are subscribed to the Google Groups "apmonitor" group.
To unsubscribe from this group and stop receiving emails from it, send an email to apmonitor+...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

--

Best regards,

John Hedengren

APMonitor Optimization Suite

http://apmonitor.com

[Ricardo LUNA HERNANDEZ]

Dear Jonh,

First of all, thanks for your time and colaboration!

Efectively, i not included the equation for xm. 

I checked your changes in the model and found an error in differential equation for x2.dt(). The term (rog*Tau2**2) that multiplicate x2.dt() is not equal or equivalent to my expression in my model code.

I included "Dcost" in my code for MVs and i could solve the problem with IPOPT :)

But i have other question;

In particual in this process, there are two constraints;

1) xm final value must be equal to 0.13 (its ok)

2) the derivate of q (quality of rice) always must be negative, i.e., q.dt() <=0 or q.dt() <0.

However, the simulations in my case and your case, q is constant in 80% and not negative over time. So, this optimal solution is not optimal for my case. It not comply this constraint, and i don't understand why not comply

Thanks in Advance.

Best regards

 

Ricardo Luna

[John Hedengren]

Ricardo,

The problem is with the parameter declaration:

K0 = m.Param(value = 1.56e+27*3600)

I recommend that you switch this to a constant in Python instead:

K0 = 3600*1.56e+27

GEKKO / APM automatically clips variables at +/- 1e20 by default, especially to encourage good model scaling. In your Intermediate expression, you use the value of K0:

K     = m.Intermediate(-K0*(x2-x1)**5*m.exp(-Ea/(R*(Tg+273.15))))

to generate K. The value of K in your prior script was very small because K0 was clipped about 1e10 below what it should have been.

K = [0.0, -1.175872e-13, 3.133979e-15, 6.470415e-14, 3.508487e-14, 5.802856e-14, 6.907941e-13, 4.253731e-12, 7.908181e-12, 8.904955e-12, 7.710384e-12, 5.867229e-12, 4.214958e-12, 2.966575e-12, 2.083597e-12, 1.472974e-12, 1.051802e-12, 7.594092e-13, 5.54308e-13, 4.087597e-13, 3.042694e-13]

The constraint q.dt()<=0 is satisfied within the tolerance that you specify for the solver. A (<) or (<=1) are equivalent for continuous variables and constraints. One of the issues why it is not decreasing is because the value of K is small. If you need q to decrease at a certain rate then I'd recommend inserting a constraint such as q.dt()<-0.

Best regards,

John Hedengren

[Ricardo LUNA HERNANDEZ]

Dear John,

Thank for your responde. I will change the parameter declaration (K0) and wil try implement contraints in q.dt()<0 or q.dt<-0.01

i have a other question about gekko;

how do you print or display the intermediates in python?

for example for case of K

because only can print variable declarations for example

var x1, x2 as;

print(x1.value)

but i cannot print intermediates.

Thank in advance!

Best regards

Ricardo

[John Hedengren]

Dear Ricardo,

 

I’m able to see the value of the intermediate (K.value). See the print statement in the following script. You may need to update your GEKKO version because we recently added the ability to print Intermediates.

 

Best regards,

 

John Hedengren

--

[Ricardo LUNA HERNANDEZ]

Dear Jonh,

Ok, i understand!

I going to update gekko. Thank a lot!

Best regards,

Ricardo