Inclure le temps en tant que variable explicite dans une contrainte dans un modèle Pyomo

Question

J'utilise PyOMO pour modéliser une réaction semi-discontinue.

Considérons un système d'ODE qui décrit un réacteur semi-discontinu où l'un des réactifs est alimenté à un débit volumique donné pour t1 unités de temps, la réaction se poursuit jusqu'à t _fin, et évidemment t1 < t _fin.

Pour spécifier l'arrêt dans le flux, je peux soit utiliser une règle conditionnelle (on suppose t1 = 3,5 * 60):

def _vol_flow_in_schedule(mod,t): 
if t<=3.5*60: 
    return mod.vol_flow_in[t] == (12.3/1000)/(3.5*60) 
else: 
    return mod.vol_flow_in[t] == 0 
m1.vol_flow_in_schedule = Constraint(m1.time,rule=_vol_flow_in_schedule) 

qui va créer une discontinuité (et mon modèle ne converge pas). Ce que je veux faire est d'utiliser une fonction sigmoïdale qui va faire passer le flux à zéro sans discontinuité.

Pour implémenter le sigmoïde, j'ai besoin de me référer à la variable de temps elle-même.

Le ci-dessous un code Matlab me donne le résultat que je veux:

t=[0:1:500]; 
acc=2; %Acceleration parameter, higher values yields sharper change. 
time_of_step=3.5*60; 
init_value = (12.3/1000)/(3.5*60); 
end_value = 0; 
sigmoidal=(init_value+(end_value-init_value)/2)... 
      +((end_value-init_value)/2)*atan((t-time_of_step)*acc)/atan(max(t)); 

Cette mise en œuvre doit cependant la variable temps explicitement dans la fonction. Comment puis-je accéder à la variable de temps à l'intérieur de la règle PyOMO? J'ai essayé le ci-dessous, mais je reçois un « ne peut pas traiter la composante scalaire « t_of_step » comme un tableau » Erreur:

m1.init_value = Param(initialize = (12.3/1000)/(3.5*60)) 
m1.end_value = Param(initialize = 0) 
m1.t_of_step = Param(initialize = 210) 
m1.acc  = Param(initialize = 5) 
. 
. 

def _vol_flow_sigmoidal (mod,t): 
return mod.vol_flow_in[t] == (mod.init_value+(mod.end_value-mod.init_value)/2)+((mod.end_value-mod.init_value)/2)*atan((t-mod.t_of_step)*mod.acc)/atan(1500) 
m1.vol_flow_sigmoidal = Constraint(m1.time,rule=_vol_flow_sigmoidal) 

Espérons que je viens de décrire clearlyt ce que je suis après. Tous les conseils sont les bienvenus,

Merci! Sal

Answer 1

Comment déclarez-vous l'index m1.time? Je suppose que vous utilisez un tableau NumPy pour initialiser l'index m1.time

Il y a un problème connu dans Pyomo (voir Issue #31) où la surcharge de l'opérateur NumPy et la surcharge de l'opérateur Pyomo finissent par se battre les uns avec les autres (NumPy se trompe en pensant que les scalaires Pyomo sont réellement indexés et tentent de les traiter comme des tableaux).

j'ai pu reproduire l'erreur avec l'exemple complet suivant:

# pyomo 4.4.1 
from pyomo.environ import *              
import numpy as np                

m1 = ConcreteModel()                
m1.time = Set(initialize=np.array([0,100,200,300,400,500]))      
m1.vol_flow_in = Var(m1.time)             

m1.init_value = Param(initialize = (12.3/1000)/(3.5*60))       
m1.end_value = Param(initialize = 0)           
m1.t_of_step = Param(initialize = 210)           
m1.acc  = Param(initialize = 5)           

def _vol_flow_sigmoidal (mod,t):             
return mod.vol_flow_in[t] == (mod.init_value+(mod.end_value-mod.init_value)/2)\ 
+((mod.end_value-mod.init_value)/2)*atan((t-mod.t_of_step)*mod.acc)/atan(1500) 
m1.vol_flow_sigmoidal = Constraint(m1.time,rule=_vol_flow_sigmoidal)    

Il y a deux alternatives qui fonctionnent, à la fois sur la base évitant d'utiliser des tableaux numpy pour initialiser ensembles Pyomo.Vous pouvez éviter soit complètement Numpy:

m1.time = Set(initialize=[0,100,200,300,400,500]) 

ou explicitement jeté le tableau NumPy à une liste:

timeArray = np.array([0,100,200,300,400,500]) 
m1.time = Set(initialize=timeArray.tolist()) 

Enfin, pour être complet, deux autres notes:

1. Ceci est également valable à initialiser ContinuousSet objets dans pyomo.dae

2. Vous verrez le même comportement même si vous évitez la déclaration explicite Pyomo Set. Autrement dit, ce qui suit va également générer l'erreur:

m1.time = np.array([0,100,200,300,400,500]) 
# ... 
m1.vol_flow_sigmoidal = Constraint(m1.time,rule=_vol_flow_sigmoidal)    

C'est parce que Pyomo va créer tranquillement l'objet Set pour vous dans les coulisses m1.vol_flow_sibmodial_index et ensuite l'utiliser ensemble pour indexer la contrainte.