Optimization problem with non-linear constraint

I am new to Optimization problems and working on a simple maximization problem, that i can solve quite simply in Excel. However, I need to scale it in Python and need some help.

I have a menu of different food items, and I need to maximize my energy output. Example:

I need to Maximize Energy(e), given the following constraints:

1. Only 1 of the food-items(i) per Macros(m) can be consumed. So I need an indicator variable (0/1) to select only 1 from each of m - Protein, Fat and Carbs.
2. Total number of calories (c) should not exceed a constant value Assume 1 portion for each item (no constraint required for this)

Problem Formulation:

Variable: X (m,i) → Binary Variable = {1 , if macro m and item i is chosen, 0 Otherwise}

Maximize e(m,i) * X(m,i)

Parameters: Calories (C) -> Calories for each (Macro, fooditem)

Subject to Constraints: For each m, Σ X (m,i) <= 1 (only 1 item per each macro can be selected) Σ c(m,i) * X(m,i)/ X(i) <= N ( calories consumer limited to constant N)

So far, I see this as a Mixed Integer Problem with a Non-linear constraint.

1. I have attempted using Pulp, but it fails due to non-linear constraint. If I remove the non-linearity, it works ok.
2. I attempted with Scipy Optimize, but Scipy doesn't allow to create integer variables.

How can I solve this using Python? Am I misinterpreting the problem here?

UPDATE: The above was missing the non-linear component that gets added due to the mean. I updated the problem from a constraint on the total to a constraint on the mean. In a non-mathematical lingo, I am taking the mean of the number I get after multiplying all the macros since I want my average calories to be less than the constant N.

So mathematically, Σ c(m,i) * X(m,i)/ X(i) <= N ( average calories consumer limited to constant N)

As you already mentioned, scipy.optimize.minimize can't handle mixed-integer problems (MIP). The most one can do is to try to solve the MIP by a penalty method, i.e. one adds a penalty function to the objective like 1.0/eps * np.sum(x*(1 - x)), where eps > 0 is a given penalty parameter and x a np.ndarray.

However, it's much more convenient to solve the problem with a MIP solver. Since your problem has a well-known knapsack-like structure, you can expect even non-commercial MIP solvers (PuLp uses CBC by default) to utilize your problem's underlying structure. Here, I'd recommend the following formulation:

Binary variables:
x[i] = 1 if fooditem i is chosen, 0 otherwise

Parameters:
a[i][m] = 1 if fooditem i covers macro m, 0 otherwise
c[i]        calories for fooditem i
e[i]        energy for fooditem i
N           total calories limit

Model:

max Σ (e[i] * a[i][m] * x[i],  ∀ i ∈ Fooditems, m ∈ Macros)

s.t. Σ (a[i][m] * x[i], ∀ i ∈ Fooditems) <= 1  ∀ m ∈ Macros. (1)
     Σ (c[i] * x[i], ∀ i ∈ Fooditems)    <= N                (2)

which can be modelled and solved like this:

import pulp

fooditems = {
    'Fish':    {'macro': 'Protein', 'calorie': 100, 'energy': 60},
    'Lamb':    {'macro': 'Protein', 'calorie': 200, 'energy': 40},
    'Egg':     {'macro': 'Protein', 'calorie': 200, 'energy': 38},
    'Banana':  {'macro': 'Carbs',   'calorie': 200, 'energy': 25},
    'Potato':  {'macro': 'Carbs',   'calorie': 200, 'energy': 30},
    'Rice':    {'macro': 'Carbs',   'calorie': 200, 'energy': 40},
    'Avocado': {'macro': 'Fat',     'calorie': 450, 'energy': 50},
    'Cheese':  {'macro': 'Fat',     'calorie': 400, 'energy': 60},
    'Cream':   {'macro': 'Fat',     'calorie': 500, 'energy': 55},
}

# parameters
macros = list({fooditems[i]['macro'] for i in fooditems})
a = {item: {m: 1 if m == fooditems[item]['macro']
            else 0 for m in macros} for item in fooditems}
c = {item: fooditems[item]['calorie'] for item in fooditems}
e = {item: fooditems[item]['energy'] for item in fooditems}
N = 1000

# pulp model
mdl = pulp.LpProblem("bla", pulp.LpMaximize)

# binary variables
x = pulp.LpVariable.dicts("x", fooditems, cat="Binary")

# objective
mdl.setObjective(sum(e[i] * a[i][m] * x[i] for m in macros for i in fooditems))

# constraints (1)
for m in macros:
    mdl.addConstraint(sum(a[i][m]*x[i] for i in fooditems) <= 1)

# constraints (2)
mdl.addConstraint(sum(x[i]*c[i] for i in fooditems) <= N)

# solve the problem
mdl.solve()

print(f"Status: {pulp.LpStatus[mdl.status]}")
for var in mdl.variables():
    print(f"{var.name} = {var.varValue:.0f}")
print(f"energy: {mdl.objective.value()}")

This yields

Status: Optimal
x_Avocado = 0.0
x_Banana = 0.0
x_Cheese = 1.0
x_Cream = 0.0
x_Egg = 0.0
x_Fish = 1.0
x_Lamb = 0.0
x_Potato = 0.0
x_Rice = 1.0
Energy: 160.0