Solve nonlinear equations without initial solution in python (matlab works, python does not)

user_11525087 0
Viewed 108

currently I try to calculate an exponential function of the following form: f(x) = aexp(-bx)+c

for context: this should be the reduction function of co2 emissions.

the conditions are following:

  • f(x=0) = 700 Mt (current co2 emissions of Germany)
  • f(x=30) = 0 Mt (no co2 emissions in 30 years)
  • integral 0->30 f(x) = budget (remaining co2 budget - area under the curve)

I can solve this problem using matlab (solve-function)

But I want to do it in python but all the solvers I tried failed. When I use the solution of the matlab script as an initial guess everything works fine. But I don't want to do that. I want to use a python solver which works if I don't know a good initial guess.

Which one should I use? I tried fsolve from scipy and gekko so far.

Here's the code for my gekko attempt:

def solve(self, end_year, budget):
        self.end_year = end_year
        self.budget = budget
        d_year = self.end_year - self.start_year
        
        m = GEKKO()
        a,b,c = [m.Var(1) for i in range(3)]
        
        eq1 = a + c - self.co2_start
        eq2 = a*m.exp(-b*d_year) + c
        eq3 = a/b*(1-m.exp(-b*d_year))+c*d_year-self.budget
        
        m.Equation([eq1==0, eq2==0, eq3==0])
        m.solve(disp=False)
        self.a, self.b, self.c = a.value, b.value, c.value
´´´
gekko scipy python
1 Answers

I haven't thought very hard about the shape of this problem, but here's a piece of code using scipy that works for me, and requires no analytical solutions or simplifications:

from math import exp
import numpy as np
import scipy.integrate
import scipy.optimize
import matplotlib.pyplot as plt

co2_start = 700
co2_end = 0
year = 30
budget = 5000 # Not mentioned in the question

def func(x, a, b, c):
    return a * exp(-b * x) + c

def func_integral(x, a, b, c):
    # In this case it's easy to work out the integral analytically, but maybe
    # you want a more sophisticated `func` later.
    return scipy.integrate.quad(func, 0, x, (a, b, c))[0]

def err_func(x):
    try:
        a, b, c = x
        err_start = func(0, a, b, c) - co2_start
        err_end = func(year, a, b, c) - co2_end
        err_area = func_integral(year, a, b, c) - budget
        return [err_start**2, err_end**2, err_area**2]
    except OverflowError:
        return [np.inf, np.inf, np.inf]

solution = scipy.optimize.fsolve(err_func, [0, 1, 0])
print(solution)

a, b, c = solution
xs = np.arange(0, year)
ys = [func(x, a, b, c) for x in xs]
plt.plot(xs, ys)
plt.show()

I'm setting a, b, c = 0, 1, 0 because we know b must be positive, and the b = 0 region is unpleasant for the solver.

It spits out an IntegrationWarning from the quad function, but I'm hoping that this is only for an extreme case far away from the actual solution.

From err_func, I'm returning an error vector here with the squared error of each of the conditions. This is a common approach, but in this case it seems to work without squaring as well.

The other ingredient is handling OverflowError because the solver tries a value for b of the order 1e6, which is too large for the exp function. We steer the solver away from this area by returning infinity.

user_146370
Related