I am trying to run a dynamic optimization with Scipy, that features a definite integral with fixed bounds. However, the optimization fails to send back any value for my objective function, and report two issues (detailed report below) :
My code is :
import numpy as np
from scipy.integrate import quad
from scipy.stats import truncpareto
from scipy.optimize import minimize, Bounds, NonlinearConstraint
φ_min = .7 # ACAP minimale
φ_max = 1.7 # ACAP maximale
p = 1 # Prix du bien manufacturier
α = 0.036 # Part de l'énergie dans la production
ν = 0.27 # Part du capital dans la production
η = 0.6 # Part du travail dans la production
ρ = 1- α - η - ν # Paramètre de rendements d'échelle décroissants
ω = 300 # Paramètre d'échelle de l'apprentissage
γ = 0.271 # Paramètre d'apprentissage
λ = 8 # Productivité verte maximal
L = 192 # Travail disponible
J = 1 # Nombre de firmes
δ = 0.1 # Taux de dépréciation
a = 2. # Paramètre de forme (shape) de la loi de Pareto
Ad = 5 # Productivité de la technologie brune
Pec = 130 # Prix de l'énergie bas carbone
Ped = 35 # Prix de l'énergie fossile
β = 0.03 # Taux d'actualisation intertemporel
CO2 = 0.4 # Facteur d'émission de l'énergie fossile
H0 = 5 # Expérience initiale verte
S0 = 100 # Budget carbone initial
T = 10 #Horizon temporel
def Integrande(phi : float ,H : float) -> float :
return (((phi*(ω**γ)+(λ*(H**γ)))/((H**γ)+(ω**γ)))**(1/ρ))*truncpareto.pdf(phi, c = φ_max/φ_min, scale = φ_min, b = a)
def EAc(H: float) -> float :
return quad(Integrande, φ_min, φ_max, args=(H))[0]
def Zt(H : float, Tax : float, Sub : float) -> float :
return ((Ad/((Ped+Tax)**α))**(1/ρ))+( EAc(H)/((1-Sub)**(ν/ρ)*Pec**(α/ρ)))
def Production(Tax : float, Sub : float, H : float) -> float:
#c = (p**(α+ν)*(α**α)*(((β*ν)/δ)**ν)*L**η*(J**(-η/ρ)))**(1/(ρ+η))*
return (Zt(H, Tax, Sub)**(ρ/(ρ+η)))*J
'transition dynamics'
def transition(States : np.ndarray , Controls : np.ndarray) -> float:
H1, S1 = States
Tax1, Sub1 = Controls
# Budget carbone
S_Next = S1 - CO2 * Zt(H1, Controls[0], Controls[1])**(-η/(ρ+η))*(Ad/((Ped + Controls[0])**(α+ρ)))**(1/ρ)*J #Fonction de conso d'énergie sale
# Expérience verte
H_Next = States[1] + (Zt(States[0], Controls[0], Controls[1])**(-η/(ρ+η)))*(EAc(Controls[1])/((1-Controls[1])**(ν/ρ)*Pec**(α/ρ)))*J
return S_Next, H_Next
'Constraint'
def contrainte(Control_seq : np.ndarray) -> float :
States = np.array([S0, H0])
Control_seq_reshaped = np.reshape(Control_seq,(T,2))
for Controls in Control_seq_reshaped:
BudgetC = transition(States,Controls)[0]
States = transition(States, Controls)
return BudgetC
'Optimization problem with objective function'
def Problème_opti(Control_seq : np.ndarray) -> float :
Prod_totale = 0
states = np.array([S0, H0])
Control_seq_reshaped = np.reshape(Control_seq,(len(Control_seq)//2,2))
t = 0
for Controls in Control_seq_reshaped:
Tax, Sub = Controls
Prod_totale = (1-β)**t * Production(Tax, Sub, states[1])
states = transition(states, Controls)
t = t+1
return -Prod_totale
def main() -> None :
result = minimize(
fun = Problème_opti,
x0 = np.zeros(2*T),
constraints = NonlinearConstraint(contrainte, lb = 0, ub = np.inf),
method = 'SLSQP'
)
return result.x, result.fun
print('optimal policy : ',main()[0])
print('Production totale', main()[1])
The console reports
c:\users\gal\dropbox\thèse\code\modèle dynamique v3.py:51: RuntimeWarning: invalid value encountered in scalar power
return ((Ad/((Ped+Tax)**α))**(1/ρ))+( EAc(H)/((1-Sub)**(ν/ρ)*Pec**(α/ρ)))
c:\users\gal\dropbox\thèse\code\modèle dynamique v3.py:66: RuntimeWarning: invalid value encountered in scalar power
H_Next = States[1] + (Zt(States[0], Controls[0], Controls[1])**(-η/(ρ+η)))*(EAc(Controls[1])/((1-Controls[1])**(ν/ρ)*Pec**(α/ρ)))*J
c:\users\gal\dropbox\thèse\code\modèle dynamique v3.py:46: IntegrationWarning: The occurrence of roundoff error is detected, which prevents
the requested tolerance from being achieved. The error may be
underestimated.
return quad(Integrande, φ_min, φ_max, args=(H))[0]
c:\users\gal\dropbox\thèse\code\modèle dynamique v3.py:64: RuntimeWarning: invalid value encountered in scalar power
S_Next = S1 - CO2 * Zt(H1, Controls[0], Controls[1])**(-η/(ρ+η))*(Ad/((Ped + Controls[0])**(α+ρ)))**(1/ρ)*J #Fonction de conso d'énergie sale
and as final value of my objective function, the program reports 'nan'
I ran the objective function ('Problème-opti') with an arbitrary argument, by calling Problème_opti(np.array([1,2,3,4,5,6,7,8,9,1])) (and commenting away the "main()" section) and obtained a rather similar feedback :
c:\users\gal\dropbox\thèse\code\untitled1.py:47: RuntimeWarning: invalid value encountered in power
return ((Ad/((Ped+Tax)**α))**(1/ρ))+( EAc(H)/((1-Sub)**(ν/ρ)*Pec**(α/ρ)))
c:\users\gal\dropbox\thèse\code\untitled1.py:63: RuntimeWarning: invalid value encountered in power
H_Next = States[1] + (Zt(States[0], Controls[0], Controls[1])**(-η/(ρ+η)))*(EAc(Controls[1])/((1-Controls[1])**(ν/ρ)*Pec**(α/ρ)))
c:\users\gal\dropbox\thèse\code\untitled1.py:42: IntegrationWarning: The occurrence of roundoff error is detected, which prevents
the requested tolerance from being achieved. The error may be
underestimated.
return quad(Integrande, φ_min, φ_max, args=(H))[0]
c:\users\gal\dropbox\thèse\code\untitled1.py:63: RuntimeWarning: divide by zero encountered in scalar divide
H_Next = States[1] + (Zt(States[0], Controls[0], Controls[1])**(-η/(ρ+η)))*(EAc(Controls[1])/((1-Controls[1])**(ν/ρ)*Pec**(α/ρ)))
I am rather confused, and I cannot see where I get this wrong. can you help me with figuring it out ?
It turns out one of the control variable should take its values between 0 and 1. The minimize function does not include such a restriction. Now I need to define a specific bound for it.