Scipy quad与optimize优化含积分和分数幂的价值函数问题求助

问题概述 #

我尝试用Scipy实现包含固定边界定积分的动态优化，但优化过程无法返回目标函数的有效值，出现两类核心问题：

• 积分舍入误差（据我理解，该误差通常在积分发散时出现）
• 标量幂运算中遇到无效值（该问题通常因对负数应用分数幂导致）

代码实现 #

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])

控制台报错信息 #

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

最终目标函数返回值为nan。

测试情况 #

调用任意参数Problème_opti(np.array([1,2,3,4,5,6,7,8,9,1]))测试目标函数，得到类似反馈：

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**(α/ρ)))

我无法定位问题所在，恳请协助排查。

内容的提问来源于stack exchange，提问作者GaL