PyOmo/Ipopt échoue avec "ne peut pas évaluer pow"

Question

J'utilise PyOmo pour générer un modèle non linéaire qui sera finalement résolu avec Ipopt. Le modèle est le suivant:

from pyomo.environ import * 
from pyomo.dae import * 

m = ConcreteModel() 

m.t = ContinuousSet(bounds=(0,100)) 

m.T = Param(default=100,mutable=True) 
m.a = Param(default=0.1) 
m.kP = Param(default=20) 
m.P = Var(m.t, bounds=(0,None)) 
m.S = Var(m.t, bounds=(0,None)) 
m.u = Var(m.t, bounds=(0,1), initialize=0.5) 

m.Pdot = DerivativeVar(m.P) 
m.Sdot = DerivativeVar(m.S) 

m.obj = Objective(expr=m.S[100],sense=maximize) 

def _Pdot(M,i): 
    if i == 0: 
    return Constraint.Skip 
    return M.Pdot[i] == (1-M.u[i])*(M.P[i]**0.75) 

def _Sdot(M,i): 
    if i == 0: 
    return Constraint.Skip 
    return M.Sdot[i] == M.u[i]*0.2*(M.P[i]**0.75) 

def _init(M): 
    yield M.P[0] == 2 
    yield M.S[0] == 0 
    yield ConstraintList.End 

m.Pdotcon   = Constraint(m.t, rule=_Pdot) 
m.Sdotcon   = Constraint(m.t, rule=_Sdot) 
m.init_conditions = ConstraintList(rule=_init) 

discretizer = TransformationFactory('dae.collocation') 
discretizer.apply_to(m,nfe=100,ncp=3,scheme='LAGRANGE-RADAU') 
discretizer.reduce_collocation_points(m,var=m.u,ncp=1,contset=m.t) 

solver = SolverFactory('ipopt') 
results = solver.solve(m,tee=False) 

Exécution des résultats du modèle dans l'erreur suivante:

Error evaluating constraint 1: can't evaluate pow'(0,0.75). 

Traceback (most recent call last): 
    File "<stdin>", line 1, in <module> 
    File "/usr/local/lib/python3.5/dist-packages/pyomo/opt/base/solvers.py", line 577, in solve 
    "Solver (%s) did not exit normally" % self.name) 
pyutilib.common._exceptions.ApplicationError: Solver (asl) did not exit normally 

La première partie de l'erreur provient de Ipopt alors que la deuxième partie vient de PyOmo. Évidemment le problème a quelque chose à faire avec le terme M.P[i]**0.75 dans les contraintes, mais changer la puissance ne résout pas le problème (bien que 2.0 ait fonctionné).

Comment puis-je résoudre ce problème?

Answer 1

Le message d'erreur indique que pow'(0,0.75) ne peut pas être évalué. Le caractère ' dans cette fonction indique la première dérivée ('' induirait la dérivée seconde). Le message dit effectivement que la première dérivée n'existe pas ou aboutit à une infinité à zéro.

Résoudre la question est facile: lié vos variables à une valeur non nulle comme suit:

m.P = Var(m.t, bounds=(1e-20,None)) 
m.S = Var(m.t, bounds=(1e-20,None))