pythonoptimizationpulp

Root mean square linearisation for linear programming

I am trying to linearise the function root mean square to use it in a linear optimisation or Mixed integer linear optimisation. Any idea how I could do this? For instance with the example below, if I wanted to maximize P*100, the model would give P=10, Q = 0 and S=10.

Many thanks

import numpy as np
import pulp

S = np.sqrt(P**2 + Q**2)

model = pulp.LpProblem("Linearise RMS", pulp.LpMaximize)

P = pulp.LpVariable("P", lowBound=-10, upBound=10 ,cat="Continuous")
Q = pulp.LpVariable("Q", lowBound=-10, upBound=10 ,cat="Continuous")
S = pulp.LpVariable("S", lowBound= 0, upBound=10 ,cat="Continuous")

objective_function  = P*100

model.setObjective(objective_function)
cbc_solver = PULP_CBC_CMD(options=['ratioGap=0.02'])
result = model.solve(solver=cbc_solver)
Solution

  • It cannot be linearised with an exact problem. It can be linearised with an approximate problem, and depending on a few things, to a decent degree of accuracy. It cannot be linearised with a continuous solver. It must be linearised with a mixed integer solver because the upper parabolic constraint is non-convex; only the lower constraint is convex. The upper constraint segments require binary selectors.

    Your toy example is poorly-chosen; the optimal solution is trivial and most of the parabolic constraints don't matter. Depending on what you're actually doing, often the lower or upper half of the parabolic constraint can and should be dropped entirely, but you've left that entirely ambiguous so I show both, despite the fact that it isn't necessary. However, I demonstrate with different bounds than your example.

    import numpy as np
import pulp
from matplotlib import pyplot as plt


def segmented_parabola(
    v: pulp.LpVariable, model: pulp.LpProblem, n_segments: int = 5, plot: bool = False,
) -> tuple[
    pulp.LpVariable,  # v**2
    list[pulp.LpVariable],  # segment selects
    plt.Axes | None,
]:
    v2 = pulp.LpVariable(name=v.name + '2', cat=v.cat)
    x = np.linspace(v.lowBound, v.upBound, n_segments)
    y = x**2
    dydx = 2*x

    if plot:
        fig, ax = plt.subplots()
        xhi = np.linspace(v.lowBound, v.upBound, 101)
        ax.plot(xhi, xhi**2)
        ax.set_xlabel(v.name)
        ax.set_ylabel(v2.name)
    else:
        ax = None

    # For the lower constraints, simply add one constraint row per line segment
    for i, (xi, yi, dydxi) in enumerate(zip(x, y, dydx)):
        model.addConstraint(
            name=f'lower_{v.name}{i}', constraint=v2 >= v * dydxi - yi,
        )
        if plot:
            ax.plot(xhi, dydxi * xhi - yi)

    # The upper constraints are non-convex, so require binary selectors
    selects = pulp.LpVariable.matrix(
        name=f'select_{v.name}', cat=pulp.LpBinary, indices=range(n_segments - 1),
    )
    model.addConstraint(
        name=f'exclusive_{v.name}', constraint=1 == pulp.lpSum(selects),
    )

    dydx = (y[1:] - y[:-1])/(x[1:] - x[:-1])
    offset = y[:-1] - dydx*x[:-1]
    # v2 <= dydxi*v + offseti + M*(1 - select)
    # M > (v2 - offseti - dydxi*v)/(1 - select)
    M = 2*(
        max(abs(v.lowBound), abs(v.upBound))**2
        - offset.min()
        - min(v.lowBound*dydx.max(), v.upBound*dydx.min())
    )

    for i, (x0, x1, dydxi, offseti, select) in enumerate(zip(
        x[:-1], x[1:], dydx, offset, selects,
    )):
        x01 = np.array([x0, x1])
        if plot:
            ax.plot(x01, dydxi*x01 + offseti)

        if x0 > v.lowBound:  # select=1 iff v >= segment left bound
            model.addConstraint(
                name=f'selectleft_{v.name}{i}',
                constraint=select <= (v - v.lowBound)/(x0 - v.lowBound),
            )

        if x1 < v.upBound:  # select=1 iff v <= segment right bound
            model.addConstraint(
                name=f'selectright_{v.name}{i}',
                constraint=select <= (v.upBound - v)/(v.upBound - x1),
            )

        # if select=1, v2 <= dydxi*v + offseti
        model.addConstraint(
            name=f'upper_{v.name}{i}', constraint=v2 <= dydxi*v + offseti + M*(1 - select),
        )

    return v2, selects, ax


def display(
    x: pulp.LpVariable, x2: pulp.LpVariable, select: list[pulp.LpVariable],
    ax: plt.Axes,
) -> None:
    print(f'{x.name} = {x.value()} ~ ±{np.sqrt(x2.value())}')
    print(f'{x2.name} = {x2.value()} ~ {x.value()**2}')
    print('selected segment', next(i for i, v in enumerate(select) if v.value()))
    print()

    ax.scatter(x.value(), x2.value())


def main() -> None:
    p = pulp.LpVariable(name='p', lowBound=-10, upBound=2.5, cat=pulp.LpContinuous)  # or -10, 10
    q = pulp.LpVariable(name='q', lowBound=-3, upBound=0.5, cat=pulp.LpContinuous)   # or -10, 10
    s = pulp.LpVariable(name='s', lowBound=0, upBound=5, cat=pulp.LpContinuous)      # or 0, 10

    model = pulp.LpProblem(name='linearise_rms', sense=pulp.LpMaximize)
    model.setObjective(p)

    p2, pa, axp = segmented_parabola(p, model, plot=True)
    q2, qa, axq = segmented_parabola(q, model, plot=True)
    s2, sa, axs = segmented_parabola(s, model, plot=True)

    model.addConstraint(name='norm', constraint=s2 == p2 + q2)
    print(model)
    model.solve()
    if model.status != pulp.LpStatusOptimal:
        raise ValueError(model.status)

    display(p, p2, pa, axp)
    display(q, q2, qa, axq)
    display(s, s2, sa, axs)
    plt.show()


if __name__ == '__main__':
    main()

    linearise_rms:
MAXIMIZE
1*p + 0.0
SUBJECT TO
lower_p0: 20 p + p2 >= -100

lower_p1: 13.75 p + p2 >= -47.265625

lower_p2: 7.5 p + p2 >= -14.0625

lower_p3: 1.25 p + p2 >= -0.390625

lower_p4: - 5 p + p2 >= -6.25

exclusive_p: select_p_0 + select_p_1 + select_p_2 + select_p_3 = 1

selectright_p0: 0.106666666667 p + select_p_0 <= 0.266666666667

upper_p0: 16.875 p + p2 + 421.875 select_p_0 <= 353.125

selectleft_p1: - 0.32 p + select_p_1 <= 3.2

selectright_p1: 0.16 p + select_p_1 <= 0.4

upper_p1: 10.625 p + p2 + 421.875 select_p_1 <= 396.09375

selectleft_p2: - 0.16 p + select_p_2 <= 1.6

selectright_p2: 0.32 p + select_p_2 <= 0.8

upper_p2: 4.375 p + p2 + 421.875 select_p_2 <= 419.53125

selectleft_p3: - 0.106666666667 p + select_p_3 <= 1.06666666667

upper_p3: - 1.875 p + p2 + 421.875 select_p_3 <= 423.4375

lower_q0: 6 q + q2 >= -9

lower_q1: 4.25 q + q2 >= -4.515625

lower_q2: 2.5 q + q2 >= -1.5625

lower_q3: 0.75 q + q2 >= -0.140625

lower_q4: - q + q2 >= -0.25

exclusive_q: select_q_0 + select_q_1 + select_q_2 + select_q_3 = 1

selectright_q0: 0.380952380952 q + select_q_0 <= 0.190476190476

upper_q0: 5.125 q + q2 + 35.875 select_q_0 <= 29.5

selectleft_q1: - 1.14285714286 q + select_q_1 <= 3.42857142857

selectright_q1: 0.571428571429 q + select_q_1 <= 0.285714285714

upper_q1: 3.375 q + q2 + 35.875 select_q_1 <= 33.21875

selectleft_q2: - 0.571428571429 q + select_q_2 <= 1.71428571429

selectright_q2: 1.14285714286 q + select_q_2 <= 0.571428571429

upper_q2: 1.625 q + q2 + 35.875 select_q_2 <= 35.40625

selectleft_q3: - 0.380952380952 q + select_q_3 <= 1.14285714286

upper_q3: - 0.125 q + q2 + 35.875 select_q_3 <= 36.0625

lower_s0: s2 >= 0

lower_s1: - 2.5 s + s2 >= -1.5625

lower_s2: - 5 s + s2 >= -6.25

lower_s3: - 7.5 s + s2 >= -14.0625

lower_s4: - 10 s + s2 >= -25

exclusive_s: select_s_0 + select_s_1 + select_s_2 + select_s_3 = 1

selectright_s0: 0.266666666667 s + select_s_0 <= 1.33333333333

upper_s0: - 1.25 s + s2 + 87.5 select_s_0 <= 87.5

selectleft_s1: - 0.8 s + select_s_1 <= 0

selectright_s1: 0.4 s + select_s_1 <= 2

upper_s1: - 3.75 s + s2 + 87.5 select_s_1 <= 84.375

selectleft_s2: - 0.4 s + select_s_2 <= 0

selectright_s2: 0.8 s + select_s_2 <= 4

upper_s2: - 6.25 s + s2 + 87.5 select_s_2 <= 78.125

selectleft_s3: - 0.266666666667 s + select_s_3 <= 0

upper_s3: - 8.75 s + s2 + 87.5 select_s_3 <= 68.75

norm: - p2 - q2 + s2 = 0

VARIABLES
-10 <= p <= 2.5 Continuous
p2 free Continuous
-3 <= q <= 0.5 Continuous
q2 free Continuous
s <= 5 Continuous
s2 free Continuous
0 <= select_p_0 <= 1 Integer
0 <= select_p_1 <= 1 Integer
0 <= select_p_2 <= 1 Integer
0 <= select_p_3 <= 1 Integer
0 <= select_q_0 <= 1 Integer
0 <= select_q_1 <= 1 Integer
0 <= select_q_2 <= 1 Integer
0 <= select_q_3 <= 1 Integer
0 <= select_s_0 <= 1 Integer
0 <= select_s_1 <= 1 Integer
0 <= select_s_2 <= 1 Integer
0 <= select_s_3 <= 1 Integer
...


Result - Optimal solution found

Objective value:                2.50000000
Enumerated nodes:               0
Total iterations:               0
Time (CPU seconds):             0.01
Time (Wallclock seconds):       0.01

Option for printingOptions changed from normal to all
Total time (CPU seconds):       0.01   (Wallclock seconds):       0.01

p = 2.5 ~ ±2.5
p2 = 6.25 ~ 6.25
selected segment 3

q = -0.375 ~ ±0.375
q2 = 0.140625 ~ 0.140625
selected segment 2

s = 2.528125 ~ ±2.5279685520195856
s2 = 6.390625 ~ 6.391416015625001
selected segment 2

    p constraints

    q constraints

    s constraints

    Depending on the size of the problem, you can "easily" (heavy quotes) scale up the resolution; this completes in less than a second with n=50:

    n=50