Problem Definition:

Pa = Constant
Pb = Constant
Vmax_a = Constant
Vmax_b = constant

Objective Function:

maximize: PaVa + PbVb - PoVo
Where: Vo = f(Po) ~ C  + e*Po


0 <= Va <= Vmax_a
0 <= Vb <= Vmax_b
Va + Vb = Vo
P_min < Po < P_max

Variables to Optimize:

Va, Vb, Po
  • 2
    $\begingroup$ What do you mean by $V_o = f(P_0) \sim C + e*P_o?$ Is this a regression equation? Is $e$ the natural base or just an arbitrary coefficient (or constant)? Is $C$ a constant? $\endgroup$
    – prubin
    Aug 1 at 15:41

2 Answers 2


Your merit function is nonlinear, continuous and differentiable.

So long as $C$ is unknown, variable and unbounded, the fact that $V_o$ may be expressed as a function of $P_o$ is immaterial. If $C$ is bounded or fixed then constraints need to be added to the demonstration below, but I assume that that is not the case.

When optimizing, $V_a$ and $V_b$ being bounded implies that $V_o$ is also bounded, but this bound can stay implicit; and the simpler form for the merit function is one where $V_o$ is substituted with $V_a + V_b$. As a cost (negative merit),

$$ c = P_o (V_a + V_b) - P_a V_a - P_b V_b $$

The gradient of the cost function is easy and can help many optimizers, so include it in your problem construction:

$$ \nabla c = \begin{bmatrix} P_o - P_a \\ P_o - P_b \\ V_a + V_b \end{bmatrix} $$

All together, an implementation could look like

import numpy as np
import pandas as pd
from scipy.optimize import check_grad, minimize

Pa = 0.3
Pb = 0.6
Pmin = 0.2
Pmax = 0.7
Vmax_a = 12
Vmax_b = 9

def cost(params: np.ndarray) -> float:
    Va, Vb, Po = params
    return Po*(Va + Vb) - Pa*Va - Pb*Vb

def cost_jacobian(params: np.ndarray) -> tuple[float, float, float]:
    Va, Vb, Po = params
    return (
        Po - Pa,
        Po - Pb,
        Va + Vb,

error = check_grad(func=cost, grad=cost_jacobian, x0=(0.9, 1.1, 1.4))
assert error < 1e-6

result = minimize(
    fun=cost, jac=cost_jacobian,
    x0=(1, 1, 1),
        (0, Vmax_a),
        (0, Vmax_b),
        (Pmin, Pmax),
assert result.success, result.message

Va, Vb, Po = result.x
Vo = Va + Vb
C = Vo - np.e*Po
df = pd.DataFrame(
        'V': (Va, Vb, Vo),
        'P': (Pa, Pb, Po),
    index=pd.Index(['a', 'b', 'o']),

print(f'\nC = {C:.3f}')
      V    P
a  12.0  0.3
b   9.0  0.6
o  21.0  0.2

C = 20.456

I am not 100% sure, but I am thinking like this: I suppose you can rewrite the objective to: $$\max V_a(P_a-P_o)+V_b(P_b-P_o)$$ since $V_a + V_b = V_o$.

And since that sum is also equal to whatever is the RHS of that strange equation @prubin asked about and since either $P_a-P_o$ or $P_b-P_o$ will always be bigger depends only on $P_a$ and $P_b$ and not on $P_o$, you can assign both $V_a$ and $V_b$ from the RHS. If $P_a > P_b$ you will max out $V_a$ before assigning the remaining to $V_b$, otherwise it's just the opposite. So you can essentially search only for $P_o$ in the range $P_\min$ and $\min(P_\max, \max(P_a, P_b))$. Especially because the RHS has to be non-negative due to being equal to the sum of two non-negative numbers -> if $P_o$ would surpass $P_a$ and $P_b$ will make the objective negative, while it may be non-negative if $P_o$ is smaller than the larger one.

The problem becomes (assuming $P_a > P_b$):

\begin{align}\max&\quad {\rm RHS} \cdot (P_a - P_o)\\ \text{s.t.}&\quad 0 \le{\rm RHS}\le {V_a}_\max + {V_b}_\max\\ &\quad P_\min \le P_o \le \min(P_\max, P_a)\end{align}

And RHS being the function of $P_o$.

Not sure how to solve that optimally though. Maybe use heuristics?


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.