Unexpected solution when using scipy.optimize.minimize to solve a optimization problem

Question

I'm trying to express the following optimization problem in scipy.

Assuming r is a known array with values: [0.96366965, 0.93341242, 0.90676733, 0.88186071, 0.8582291, 0.83472442, 0.80977363, 0.783683, 0.75743122, 0.73259614]

def objective_function(x):
    e_t = np.zeros(10)
    e_t[0] = 1 - 1 / (1 + r[0])

    for t in range(1, 9):
        sub = 0
        for ti in range(0, t):
            sub += np.sum(x[: (ti + 1)]) / (1 + r[ti]) ** (ti + 1)

        e_t[t] = 1 - 1 / (1 + r[t]) ** (t + 1) - sub

    e_t[9] = 0

    return np.sum(e_t**2)

def constraint_function(t):
    def calculate_et(x):
        sub = 0
        for ti in range(0, t):
            sub += np.sum(x[: (ti + 1)]) / (1 + r[ti]) ** (ti + 1)

        e_t = 1 - 1 / (1 + r[t]) ** (t + 1) - sub
        return e_t

    return calculate_et

x0 = np.zeros(10)
cons = (
    {"type": "ineq", "fun": constraint_function(1)},
    {"type": "ineq", "fun": constraint_function(2)},
    {"type": "ineq", "fun": constraint_function(3)},
    {"type": "ineq", "fun": constraint_function(4)},
    {"type": "ineq", "fun": constraint_function(5)},
    {"type": "ineq", "fun": constraint_function(6)},
    {"type": "ineq", "fun": constraint_function(7)},
    {"type": "ineq", "fun": constraint_function(8)},
    {"type": "ineq", "fun": constraint_function(9)},
)

result = minimize(
    objective_function,
    x0,
    bounds=Bounds(lb=0),
    constraints=cons,
)
alpha = result.x

c = np.zeros(10)
c[0] = alpha[0]
for i in range(1, 10):
   c[i] = alpha[i-1] + alpha[i]

However, the solution I got was quite off from the expected values, so I'm assuming this is not the right method to solve the problem. What is wrong with my current approach and are there other ways to approach it?

the solution I got for c is [9.02212140e-01 9.02212140e-01 2.39808173e-13 1.06359366e-13 4.57411886e-14 1.75415238e-14 6.16173779e-15 2.17534324e-15 4.67095103e-16 1.25975522e-17]

the expected solution is closer to [0.0318, 0.0318, 0.0318, 0.0318, 0.0318, 0.0318, 0.0320, 0.0325, 0.0336, 0.0353]

Answer 1

Check your math. With your current objective, the "expected solution" has a much higher objective cost than what minimize is able to come up with:

from functools import partial

import numpy as np
from scipy.optimize import minimize, Bounds, NonlinearConstraint


def calculate_e(alpha: np.ndarray, t: np.ndarray, r: np.ndarray) -> np.ndarray:
    den = 1/(1 + r[:-1])**t
    inner = alpha[:-2].cumsum() * den[:-1]
    e19 = 1 - den - np.concatenate(([0], inner.cumsum()))
    return np.concatenate((e19, [0]))


def objective_function(alpha: np.ndarray, t: np.ndarray, r: np.ndarray) -> float:
    e = calculate_e(alpha, t, r)
    return e.dot(e)


def main() -> None:
    r = np.array((
        0.96366965, 0.93341242, 0.90676733, 0.88186071, 0.85822910,
        0.83472442, 0.80977363, 0.78368300, 0.75743122, 0.73259614,
    ))
    n = r.size
    t = np.arange(1, n)
    c_expected = np.array((
        0.0318, 0.0318, 0.0318, 0.0318, 0.0318,
        0.0318, 0.0320, 0.0325, 0.0336, 0.0353,
    ))
    alpha_expected = np.concatenate((
        c_expected[:1],
        np.diff(c_expected),
    ))

    result = minimize(
        fun=objective_function,
        x0=alpha_expected,
        args=(t, r),
        bounds=Bounds(lb=0),
        constraints=NonlinearConstraint(
            fun=partial(calculate_e, t=t, r=r),
            lb=0, ub=np.inf,
        ),
        tol=1e-9,
    )
    assert result.success, result.message
    alpha = result.x
    c = alpha.cumsum()

    print(f' Initial cost: {objective_function(alpha_expected, t, r):.2f}')
    print(f'Solution cost: {result.fun:.2f}')
    print('α =')
    print(alpha)
    print('C =')
    print(c)
    print('e =')
    print(calculate_e(alpha, t, r))


if __name__ == '__main__':
    main()

 Initial cost: 6.71
Solution cost: 0.35
α =
[9.05392583e-01 1.38775468e-16 4.16322395e-17 2.77550652e-17
 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00
 1.10000000e-03 1.70000000e-03]
C =
[0.90539258 0.90539258 0.90539258 0.90539258 0.90539258 0.90539258
 0.90539258 0.90539258 0.90649258 0.90819258]
e =
[4.90749373e-01 2.71411502e-01 1.52473537e-01 8.63851746e-02
 4.87947275e-02 2.68482231e-02 1.36020676e-02 5.32970702e-03
 2.22044605e-16 0.00000000e+00]