Question: 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) 

Answer: 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     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: