Can scipy.optimize Find Optimal Input Values When Multiple Products are Involved?

I'd like to find the optimal values for Input A for Product 1 and for Input A for Product 2 with the aim to maximize Total Output and subject to a given constraint. I've tried using Python's Scipy minimize function and it works if I have just one Product but it does not work for multiple Products.

Can scipy.optimize Find Optimal Input Values When Multiple Products are Involved?

Here is what I tried assuming just two Products exist (In reality I have several thousand such Products):

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

Product1_InputB = np.array([0.5])
Product1_InputC = np.array([1])
Product1_InputD = np.array([1])
Product1_InputE = np.array([0.08])
Product1_InputF = np.array([20])

Product2_InputB = np.array([0.5])
Product2_InputC = np.array([1])
Product2_InputD = np.array([2])
Product2_InputE = np.array([0.1])
Product2_InputF = np.array([30])

def Neg_Product1_Output(Product1_InputA):
    return -1 * ((2.71828**((Product1_InputA-Product1_InputB)*(0.5*Product1_InputC-1.5*Product1_InputD)))/(1+(2.71828**((Product1_InputA-Product1_InputB)*(0.5*Product1_InputC-1.5*Product1_InputD))))*(Product1_InputA-Product1_InputB))

def Neg_Product2_Output(Product2_InputA):
    return -1 * ((2.71828**((Product2_InputA-Product2_InputB)*(0.5*Product2_InputC-1.5*Product2_InputD)))/(1+(2.71828**((Product2_InputA-Product2_InputB)*(0.5*Product2_InputC-1.5*Product2_InputD))))*(Product2_InputA-Product2_InputB))

def Neg_Total_Output(Product1_InputA,Product2_InputA):
    return Neg_Product1_Output + Neg_Product2_Output

def constraint(Product1_InputA, Product2_InputA):
    return (((Product1_InputA - Product1_InputE) * Neg_Product1_Output) + ((Product2_InputA - Product2_InputE) * Neg_Product2_Output)) / Neg_Total_Output - 2

con = {'type':'ineq', 'fun': constraint}

Product1_InputA_Initial_Guess = np.array([3])
Product1_InputA_Initial_Guess = np.asarray([3])
Product2_InputA_Initial_Guess = np.array([1])
Product2_InputA_Initial_Guess = np.asarray([1])
Product1_bound = [(0.3,4)]
Product2_bound = [(0.3,4)]

optimized_results = minimize(Neg_Total_Output,Product1_InputA_Initial_Guess,bounds=Product1_bound,constraints=con)

Product1_InputA_Optimal = optimized_results.x
Product1_InputA_Optimal

When I run the line optimized_results = ... I get the below error:

TypeError: constraint() missing 1 required positional argument: 'Product2_InputA'

Am not sure how to include Product2 in the optimized_results minimize function above.

You have stability and convergence issues, but at least to address the "multiple products" question this is very straightforward:

from functools import partial

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


def neg_product_output(
    a: np.ndarray,
    b: np.ndarray,
    c: np.ndarray,
    d: np.ndarray,
    e: np.ndarray,
    f: np.ndarray,
) -> np.ndarray:
    inner = np.exp(
        (a - b)*(0.5*c - 1.5*d)
    )
    return -inner/(1 + inner)*(a - b)


def neg_total_output(
    a: np.ndarray,
    bcdef: np.ndarray,
) -> float:
    out = neg_product_output(a, *bcdef)
    return out.sum()


def constraint(
    a: np.ndarray,
    bcdef: np.ndarray,
) -> float:
    out = neg_product_output(a, *bcdef)
    b, c, d, e, f = bcdef
    return (a - e).dot(out) / out.sum()


def solve(
    bcdef: np.ndarray,
) -> np.ndarray:
    # con - 2 >= 0: con >= 2
    con = NonlinearConstraint(
        fun=partial(constraint, bcdef=bcdef),
        lb=2,
        # The solver does not support an infinite upper bound to this constraint
        ub=100,
    )

    a12_guess = np.full(shape=bcdef.shape[1], fill_value=3.)
    lb = np.full_like(a12_guess, fill_value=0.3)
    ub = np.full_like(a12_guess, fill_value=2.25)
    f = bcdef[-1, :]
    ub[f == 30] = 2.5
    ub[f >= 40] = 3.0

    results = minimize(
        fun=neg_total_output, args=(bcdef,),
        x0=a12_guess, bounds=Bounds(lb=lb, ub=ub),
        constraints=con, tol=1e-12,
    )
    if not results.success:
        raise ValueError(results.message)
    print(results.message, 'in', results.nit, 'iterations')
    return results.x


def main() -> None:
    a = solve(
        bcdef=np.array((
            ( 0.5 ,  0.5 ),
            ( 1.  ,  1.  ),
            ( 1.  ,  2.  ),
            ( 0.08,  0.1 ),
            (20.  , 30.  ),
        )),
    )
    print('a =')
    print(a)


if __name__ == '__main__':
    main()

Optimization terminated successfully in 12 iterations
a =
[2.25       1.51094911]