maximize (with scipy) production having constrained budget

Question

I tried to solve common optimization problem: to maximize profit from certain production function having limited budget (7800$), but still something is wrong in such logics of optimization:

from scipy.optimize import minimize

def max_out(x):     # Q_func
    return 36*x[0]**(1/2) * x[1]**(1/3) * x[2]**(1/4)

def obj(x):     # maximize total P
    y= -(25*x[0] + 20*x[1] + 10*x[2])  # maximize production revenue
    return y

def constr(x):
    return ( 7800- (25*x[0] + 20*x[1] + 10*x[2])  )  # constraint is budget 7800

cons = ({'type': 'ineq', 'fun':  constr },
        {'type': 'ineq', 'fun': lambda x: x[0] },
        {'type': 'ineq', 'fun': lambda x: x[1] },
        {'type': 'ineq', 'fun': lambda x: x[2] })

##bnds = ((0, None), (0, None), (0, None))      # bounds= bnds,
res = minimize(obj, (10, 10, 10), method='SLSQP',  constraints=cons)
print(res.x)
r= [round(x) for x in res.x]

print("raw materials needed for production func:", r)
print(f'to maximize production revenue to {-obj(r):.2f} $')
print("cost of product: ", max_out(r))

# raw materials needed for production func: [171, 139, 74]
# to maximize production revenue to 7795.00 $
# cost of product:  7152.316966012599

as is resulted: I've, probably, got the closest to budget total_price of produced goods from given raw materials [x0, x1, x2]... but max_out func gives cheaper total_price for all produced. But I need max_out being close to the budget (as is production price), in order to find the price for sale (overall cost, that should be higher than inputted budget)... something is wrong!.. or have I formulated the task to the python in an incorrect way ??

p.s. frankly speaking I'm getting smaller price of product compared to the raw materials not for the first time, while trying to solve tasks like this - but this sounds strange for me... what is my misunderstanding ? and how to change the code to input raw_materials to budget & maximize the total_cost of produced units ?

Answer 1

for speed for huge data I suppose stochastic method as well (with this advice) suits, and even without revenue func (as final costs are inputed in bounds) - no func to calc further res.fun gives Y relative used bounds. In general I think there is a Dual Problem (minimizing production & maximizing product cost - thus can be solved either: either Maximazing output subject to budget, or Minimize cost subject to budget, hence define objective func correspondingly).

And besides

The differential evolution method is stochastic in nature. It does not use gradient methods to find the minimum, and can search large areas of candidate space, but often requires larger numbers of function evaluations than conventional gradient-based techniques.

import numpy as np
from scipy.optimize import  differential_evolution
from scipy.optimize import  LinearConstraint

unit_costs = [25,20,10]
total_budget = 78_000

def max_out(x):     # Q_func
    return -(36*x[0]**(1/2) * x[1]**(1/3) * x[2]**(1/4))

cons = LinearConstraint( unit_costs, lb=0, ub=total_budget)

bnds = ((0, total_budget/unit_costs[0]),
        (0, total_budget/unit_costs[1]),
        (0, total_budget/unit_costs[2]))

result = differential_evolution(max_out,  bounds=bnds, constraints=cons,  seed=1, tol= 1)

print("raw materials limited to budget:", list(map(lambda x: x.round(), result.x)))
print("revenue of production func maximized to budget:", -max_out(result.x).round())
print("revenue of production func maximized to budget:", -result.fun.round())

though stochasticity of algo can give only probabilistic guarantee of recieving the result, I suppose, but adds the speed & is more suitable for finding Global minimum (as well as simulated annealing)... perhaps, it even lacks the problem of stucking in local_minimum ?? (I'm not sure?).

but usually DE is An alternative to nonlinear convex optimization , as stated at link

P.S. I just don't understand from this description - What is the difference between Differential_Evolution & Genetic_Algo ??

Answer 2

Don't make minimize infer a numerical Jacobian when you can provide one yourself.

Unless you have a really good reason, you can leave the method as the default.

If your raw material quantities need to be integral, then minimize is not the right approach for this problem and you need MINLP. In this case you got lucky because the optimal solution happens to have numbers that are nearly integral anyway. In fact, with the other adjustments in this answer the accuracy of the output is improved to four decimal digits. (You could narrow tol if you cared to improve this further, though that probably isn't needed.)

For the orders of magnitude in your question, (1, 1, 1) is not a good initial guess, and (1000, 1000, 1000) is more appropriate.

Rather than needing to double-negate your budget, split the positive and negative-sense budget expressions into separate functions.

import numpy as np
from scipy.optimize import check_grad, minimize, LinearConstraint

powers = 1/np.arange(2, 5)
revenue_c = np.array((25, 20, 10))
max_budget = 78_000


def budget(x: np.ndarray) -> float:
    return 36 * (x**powers).prod()


def maximize_budget(x: np.ndarray) -> float:
    return -budget(x)


def jacobian(x: np.ndarray) -> np.ndarray:
    x0, x1, x2 = x
    p0, p1, p2 = x**powers
    return -36 * powers * (
        x0**-0.5   * p1*p2,
        x1**(-2/3) * p0*p2,
        x2**-0.75  * p0*p1,
    )


def max_revenue(x: np.ndarray) -> float:
    return revenue_c.dot(x)


x0 = (1_000,)*3
error = check_grad(func=maximize_budget, grad=jacobian, x0=x0)
assert error < 1e-3

res = minimize(
    fun=maximize_budget, jac=jacobian, x0=x0,
    bounds=((0, max_budget),)*3,
    constraints=LinearConstraint(A=revenue_c, ub=max_budget),
)
assert res.success, res.message

np.set_printoptions(precision=2)
print('raw materials needed for production func:', res.x)
print(f'to input cost of raw materials into budget ${max_revenue(res.x):,.2f}')
print(f'cost of product: ${budget(res.x):,.2f}')

'''
raw materials needed for production func: [1440. 1200. 1800.]
to input cost of raw materials into budget $78,000.00
cost of product: $94,557.42
'''

Answer 3

Adding to your solution, I advise you use a LinearConstraint object since your constraint is linear, this gives more accurate computation of its derivatives. Plus defining bounds is usually a good idea as it helps greatly the solver in carving out the space in which to search on. See my proposed solution (heavily copied on yours):

from scipy.optimize import minimize, LinearConstraint
import numpy as np

unit_costs = [25,20,10]
total_budget = 78_000

def max_out(x):     # Q_func
    return -(36*x[0]**(1/2) * x[1]**(1/3) * x[2]**(1/4))

def max_revenue(x):     # maximize total P
    return -np.dot(x, unit_costs)


cons = LinearConstraint( unit_costs, lb=0, ub=total_budget)

bnds = ((0, total_budget/unit_costs[0]),
        (0, total_budget/unit_costs[1]),
        (0, total_budget/unit_costs[2]))
res = minimize(max_out, x0=(1,1,1), method='SLSQP',  bounds=bnds, constraints=cons)
assert res.success, res.message
print(res.x)
r= [round(x) for x in res.x]

print("raw materials needed for production func:", r)
print(f'to input cost of raw materials into budget {-max_revenue(r):.2f} $')
print("cost of product: ", -max_out(r))

Answer 4

I finally got correct results (compared to manual calculus - p9 here):

from scipy.optimize import minimize

def max_out(x):     # Q_func
    return -(36*x[0]**(1/2) * x[1]**(1/3) * x[2]**(1/4))    # maximize to budget limit

def max_revenue(x):     #  total P
    y= (25*x[0] + 20*x[1] + 10*x[2])  #  product P
    return y

def constr(x):
    return ( 78000- (25*x[0] + 20*x[1] + 10*x[2])  )  # constraint is budget 7800

cons = ({'type': 'ineq', 'fun':  constr },
        {'type': 'ineq', 'fun': lambda x: x[0] },
        {'type': 'ineq', 'fun': lambda x: x[1] },
        {'type': 'ineq', 'fun': lambda x: x[2] })

##bnds = ((0, None), (0, None), (0, None))      # bounds= bnds,
res = minimize(max_out, (10, 10, 10), method='SLSQP',  constraints=cons)
print(res.x)
r= [round(x) for x in res.x]

print("raw materials needed for production func:", r)
print(f'to input cost of raw materials into budget {max_revenue(r):.2f} $')
print("cost of product: ", -max_out(r))

I was mistaken in used func for optimization - correct is: res = minimize(max_out, (10, 10, 10), method='SLSQP', constraints=cons)