Multi-objective LP with PuLP in Python

In some of my posts I used lpSolve or FuzzyLP in R for solving linear optimization problems. I have also used PuLP and SciPy.optimize in Python for solving such problems. In all those cases the problem had only one objective function.

In this post I want to provide a coding example in Python, using the PuLP module for solving a multi-objective linear optimization problem.

A multi-objective linear optimization problem is a linear optimization problem with more than just one objective function. This area of linear programming is also referred to as multi-objective linear programming or multi-goal linear programming.

Below I stated an examplaric multi-objective linear optimization problem with two objective functions:

Assuming that in above problem statement the two objective functions represent two different goals, such as e.g. service level and profit margin of some product portfolio, I test two alternative approaches for solving this problem.

The first approach will be to solve for one of the objectives, then fix the problem at the optimal outcome of that first problem by adding an additional constraint to a second optimization run where I will then maximize the second objective function (subject to the constraint of keeping the optimal objective value to the first sub-problem).

The second approach will be to add the two objectives together, i.e. to merge them into one objective function by applying weights. By sampling the weights and solving the combined problem for each sampled weight the optimal outcome can be reviewed in dependency of the weights.

Approach 1: Maximizing for one objective, then adding it as a constraint and solving for the other objective

Using PuLP I maximize the first objective first, then I add that objective function as a constraint to the original problem and maximize the second objective subject to all constraints, including that additional constraint.

In mathematical syntax, the problem we solve first can be stated as follows:

Here is the implementation of above problem statement in Python, using the PuLP module:

# first, import PuLP
import pulp

# then, conduct initial declaration of problem
linearProblem = pulp.LpProblem("Maximizing for first objective",pulp.LpMaximize)

# delcare optimization variables, using PuLP
x1 = pulp.LpVariable("x1",lowBound = 0) 
x2 = pulp.LpVariable("x2",lowBound = 0) 

# add (first) objective function to the linear problem statement
linearProblem += 2*x1 + 3*x2 

# add the constraints to the problem
linearProblem += x1 + x2 <= 10
linearProblem += 2*x1 + x2 <= 15

# solve with default solver, maximizing the first objective
solution = linearProblem.solve()

# output information if optimum was found, what the maximal objective value is and what the optimal point is
print(str(pulp.LpStatus[solution])+" ; max value = "+str(pulp.value(linearProblem.objective))+
      " ; x1_opt = "+str(pulp.value(x1))+
      " ; x2_opt = "+str(pulp.value(x2)))
Optimal ; max value = 30.0 ; x1_opt = 0.0 ; x2_opt = 10.0

Now, I re-state the original problem such that the second objective function is maximized subject to an additional constraint. That additional constraint requests that the first objective must be at least 30. Using mathematical syntax the problem I now solve can be stated as follows:

Here is the implementation of above problem statement in Python, using PuLP:

# remodel the problem statement
linearProblem = pulp.LpProblem("Maximize second objective",pulp.LpMaximize)
linearProblem += 4*x1 - 2*x2
linearProblem += x1 + x2 <= 10
linearProblem += 2*x1 + x2 <= 15
linearProblem += 2*x1 + 3*x2 >= 30

# review problem statement after remodelling
linearProblem
Maximize_second_objective:
MAXIMIZE
4*x1 + -2*x2 + 0
SUBJECT TO
_C1: x1 + x2 <= 10

_C2: 2 x1 + x2 <= 15

_C3: 2 x1 + 3 x2 >= 30

VARIABLES
x1 Continuous
x2 Continuous

Now, I solve this problem, using the default solver in PuLP:

# apply default solver
solution = linearProblem.solve()

# output a string summarizing whether optimum was found, and if so what the optimal solution 
print(str(pulp.LpStatus[solution])+" ; max value = "+str(pulp.value(linearProblem.objective))+
      " ; x1_opt = "+str(pulp.value(x1))+
      " ; x2_opt = "+str(pulp.value(x2)))
Optimal ; max value = -19.999999999995993 ; x1_opt = 1.0018653e-12 ; x2_opt = 10.0

This approach suggests that x1 = 0 and x2 = 10 is the the optimal solution. The optimal objective values would be 30 for objective one, and -20 for objective two.

Approach 2: Combine objectives, using sampled weights and iterations with defined step-size

When applying this approach, we will restate the original problem as follows:

The question now is how to choose α.

A typical approach in a situation like this is to identify an efficient frontier. In economics this is e.g. known as “pareto optimality”. To construct such an approach I sample alpha in steps of 0.01. For each value of alpha I restate the problem, using PuLP – and solve it.

I store my results in a list and visualize the outcome using matplotlib.pyplot:

# import matplotlib.pyplot
import matplotlib.pyplot as plt 

# import pandas and numpy for being able to store data in DataFrame format
import numpy as np
import pandas as pd

# define step-size
stepSize = 0.01

# initialize empty DataFrame for storing optimization outcomes
solutionTable = pd.DataFrame(columns=["alpha","x1_opt","x2_opt","obj_value"])

# iterate through alpha values from 0 to 1 with stepSize, and write PuLP solutions into solutionTable
for i in range(0,101,int(stepSize*100)):
        # declare the problem again
        linearProblem = pulp.LpProblem("Multi-objective linear maximization",pulp.LpMaximize)
        # add the objective function at sampled alpha
        linearProblem += (i/100)*(2*x1+3*x2)+(1-i/100)*(4*x1-2*x2)
        # add the constraints
        linearProblem += x1 + x2 <= 10
        linearProblem += 2*x1 + x2 <= 15
        # solve the problem 
        solution = linearProblem.solve()
        # write solutions into DataFrame
        solutionTable.loc[int(i/(stepSize*100))] = [i/100,
                                                     pulp.value(x1),
                                                     pulp.value(x2),
                                                     pulp.value(linearProblem.objective)]

# visualize optimization outcome, using matplotlib.pyplot
# -- set figure size
plt.figure(figsize=(20,10))
# -- create line plot
plt.plot(solutionTable["alpha"],solutionTable["obj_value"],color="red")
# -- add axis labels
plt.xlabel("alpha",size=20)
plt.ylabel("obj_value",size=20)
# -- add plot title
plt.title("Optimal combined objective function value as a function of alpha",size=32)
# -- show plot
plt.show()

To complete this article I print out the head of the optimization outcome DataFrame table:

solutionTable.head()
alphax1_optx2_optobj_value
00.007.50.030.00
10.017.50.029.85
20.027.50.029.70
30.037.50.029.55
40.047.50.029.40

Related content

If you are interested in mathematical programming and optimization in Python you can consider the following courses and programs:

You May Also Like

Leave a Reply

54 comments

ahmed says:

hello, i want to know how can i get all the solutions of linear optimization problem using pip ?

Achyu says:

What if we have one objective to maximize and the other objective to minimize? How would that change the logic used in the article?

Hi Achyu,
thank you for your question.
I would change the direction of the minimization objective. So if e.g. the objective to minimize is x1+x2, then change the direction to -x1-x2 and maximize it. In that way you could apply the same logic as in this article.

If you post an exemplary problem here, I will model and solve it for you.

Junnu says:

Hi Linnart, can you please solve this
Consider the following biobjective integer knapsack problem (BOIKP):
maximize f(x) = [f1(x) = 8×1 + 9×2 + 3×3, f2(x) = 3×1 + 2×2 + 10×3]
subject to 2×1 + 2×2 + 3×3 ≤ 6
x1, x2, x3 ≥ 0, integer

Hi, If you need help with concrete model implementations or problem implementations you can contact us via the contact us form and hire a consultant for a coaching session (paid service). Thank you and have a nice day.

Fatima Tahir says:

Thanks a lot it was so helpful

Kishalay Bairagi says:

instead of linearProblem += 2*x1 + 3*x2 can a function be defined as objective function like linearProblem += F(x1, x2) where (2*x1 + 3*x2) is being calculated inside the function F(x1, x2)

Dear Kishalay. Excellent question. It is possible to do this, yes! Referring to your example, write the function like this:

def f1(x1,x2):
return 2*x1 + 3*x2

and then write

linearProblem += f1(x1,x2)

Kishalay Bairagi says:

I have one more question. The running time of MOILP which produces exact method is far less than the running time of pareto based MOILP method. Therefore, while comparing the running time of MOILP based approach with a meta-heuristic based approach (like NSGA-II) should I use the running time of exact or pareto based MOILP method???

Kishalay, thank you for your question. What do you mean with “pareto-based” and what do you mean with “exact”. It is not clear to me. Thank you for elaborating further.

Anjani Keshari says:

Got a problem….Linnart
here it is …

Both objectives have scaling issues …

1st Obj – Sum of quantities supplied by supplier (Maximization obj)

2nd Obj – quantities provided by Supplier * Cost quoted by them (Min Obj)

here second objective has cost and that is causing scaling issue and overall objective becomes negative (Maximize ) and hence overall becomes minimization problem

thus i need to use probability for first objective >0.95 in order to get output else output is coming as zero

Is there a simplified way to handle scaling of blended objective (multi objective problems)

Hello Anjani.

Thank you for your question.

In this post I demonstrated two possible approaches:
1) adding the first objective as a constraint and optimize the second objective, with the constraint of achieving the first objective to a specified extend
2) adding objectives, resulting in a joint weighted objective function

In this case I would use the first approach. That is, solve e.g. the problem of minimizing the costs while adding the constraint of sourcing some minimum amount bz the respective supplier. Depending on the complexity you will have to set up many sub-problems, with different lower values for your “sourcing” constraints.

And with regards to maximization vs. minimization: If you minimize costs, it is the same as maximizing “negative costs”. For example instead of solving min x1+x2, solve max -x1-x2 instead.

Anjani Keshari says:

Hi Linnart,

Very thanks for the reply.
But i need weights to play with trade off between both objectives which is not their in first method.

I was going through a paper. Kindly see section 2.1.1. It talks about taking (max-min) value of each objective to divide each objective and then use weights.

Do you think this may help in scaling objectives in order to say when alpha = 0.5 and (1-alpha) = 0.5 then its equal weightage for both the objectives ?

Hi Anjani.
Thanks for the questions and discussion. A couple of thoughts on your comments:
1) The first approach can in my opinion be used to play with trade-offs. Because you can set the contraints values differently, depending on how much you are willing to “sacrifice” the given objective. It off course means you would have to solve the problem many times, with different constraint values for the objective added as a constraint. It is a bit of work and a little messy, but I still think it can solve your issue. The idea is essentially the same as the one you cite from section 2.1.1 – which I comment below in 2).
2) The idea of “normalizing” the objective function by dividing it through min and max essentially means that when you weigh them with a coefficient alpha in the range 0% to 100% you are assigning priorities to the unconstrainted optimum of each single objective. It is an effective way of normalizing in this case and will allow you to investigate the sensitivity of both objectives with regards to the solution space.

I am not sure whether the approach discussed in my above comment 1) will give you better insights at less effort than 2), or vice versa. It would be interesting to hear about your findings. I will try to make an example that addresses your question when I find the time to do so.

Thank you for sharing this paper. As I mentioned in my other comment I believe this is effective and a good approach. But the discussed approach #1 in this post, adding one of the objectives as a constraint, can in my opinion also be used for playing with trade-offs – by solving the problem over and over again at different constraint values.

What I am not sure of is which approachs is fastest and most effective in terms of providing insights – for your specific problem and for all linear problems in general. I am not sure about whether one of the approaches might be better in some situations but “worse” (i.e. inferior) in other situations.

Anjani says:

Hi Linnart,

As discussed, i have few observations from both the approaches i tried.

Approach1 –
I tried to multiply alpha (from 0 till 1) with first objective and used it as constraint with second objective. As you told it suggests how much one wants to sacrifice. Problem here is user do not know what and how much it wants to sacrifice. Also since for a set of (w1, w2), it will have to solve two objectives one by one, it tends to take more time than Approach 2.

Approach 2-
It suits our problem. We see blending considers multi objective as single objective, which is faster than App – 1 and plays b/w the trade off b/w two objectives. Only one issue i found that in one of the test cases, even after doing normalization for each objective with (max-min), we saw that at weights (0.5, 0.5), though it should consider both objective as equally important, first objective was more prevalent than than second. Is there a way to handle it ?

Hi again Anjani, just to understand you clearly – when you say normalizing the objective with max-min, what do you mean exactly? You normalize the objective by dividing through (a) the maximum or through (b) the maximum less minimum? I assume (a), right?

Anjani Keshari says:

Hi Linnart, Sorry for the confusion. Its only Max not (Max – min). Thanks for correcting.

Hi Anjani, just because alpha is equal to 50% does not mean the relative objective achievement of both objectives has to be the same. And also the solution will not necessarily be the “middle” of the two isolated objectives. For exampe f1(x1,x2) = 2×1 + 2×2 is objective 1; f2(x1,x2) = 3×1-4×2 is objective 2. Constraint it x1,x2 <= 10 and x1,x2>=0. Optimal solution to f1 is x1=10,x2=10 with f1=40; optimal solution to f2 is x1=10,x2=0 with f2 = 30. If we combine the normalized objective functions f1/40 and f2/30 with alpha=50% then the optimal solution is x1=10,x2=0. In that case f1 is achieved by 50%, f2 is achieved by 100%.

Anjani Keshari says:

But then the problem arises that how to explain the choice of weightage to the client/user. It’s confusing and is there a way to handle it. Or there are alternative methods as it seems to be a flaw to me right now.

Hi Anjani, I think what you are asking for is to achieve both objectives by 50% relative to their isolated optimum, correct? How do you know that this is possible? It might not be possible to find a solution like that. This is also why I like the other approach more: Solve for the first objective, then add it as a constraint to a second version of the problem focusing on the second objective. E.g. you could say objective function value no. 1 should achieve a value equal to 50% of the isolated optimum.

andrei nunes says:

Hello Linnart!!
Thanks for share this approach in pulp, im trying to model and solve a personel schedduling / rostering problem right now. I will use the firts approach, becouse a i want to minimize cost and maximize the satisfaction of the workwers.
My highly priority is to minimise the cost, so i think this will work well. My problem is NP-Hard, so it could be difficult to solve to a large “n”.
Do you know some euristhica to reduce the time of the solution in pulp?
Andrei Nunes

I would maybe need some more information on what problem it is and how you implemented it. But in general, for very hard problems it can be interesting to apply local search optimization combined with some greedy-based improvement strategy.

Andrei Nunes says:

Thanks for your response!

I will study a little more about it.

Regards!!

You are welcome. If you have a more specific problem description then you are always welcome to post it here in the comment section 🙂

Anjani Kumar Keshari says:

Hi Linnart,

Thanks for all the suggestions you made. It worked with our work and we were able to complete our first phase.

I am now working on another project as well ‘Vehicle Route Problem’. It’s an interesting one.

You are a deserving Guru. Will keep looking such support and guidance from you.

Thanks,
Anjani

Hi Anjani. For vehicle routing problem maybe also check Google OR-tools documentation, e.g. for Python. They got some nice vehicle routing problem examples and demonstrate how to solve it with Google OR-tools

Anjani Kumar Keshari says:

Hi Linnart,

Yes for now we are using Google Or Tools. But i am afraid it doesn’t handle multiple Objective. As in we have to minimize distance as well as maximize delivery of items from site visits. Is there any way to handle with Google OR or other approach.

Thanks,
Anjani

Sai says:

Hi Linart,
Thank you for sharing this approach. I have one more question. How do you solve if there are more than 2 equations to solve. alpha and 1-alpha grid will increase exponentially. How do you think we can tackle this one?

Anjani Kumar Keshari says:

Hi Linnart,

Yes for now we are using Google Or Tools. But i am afraid it doesn’t handle multiple Objective. As in we have to minimize distance as well as maximize delivery of items from site visits. Is there any way to handle with Google OR or other approach.

Thanks,
Anjani

Anjani Kumar Keshari says:

Sorry mistakenly replied to your comment.

Hello Sai, thank you for your question. Sorry for only seeing it now, as I thought your question had been answered by someone else.

For multi-objective optimization problems there is not a “single-true-approach”. Meaning, due to multi-objective optimization being closely related to the concept of pareto-optimality, there is not a standard procedure for how to solve problems with multiple conflicting objectives.

If objectives are not in conflict, it becomes easier and once could more easily scalarize them into a single objective. But usually such problems have conflicting objectives.

The approaches defined in this article can be used for more than two objectives as well. When e.g. scalalarizing into a single objective using weight-factors, once can add additional parameters, like beta, gamme etc. and then create a table overview of optimal solutions to the scalarized problem subject to defined weight values. A graphical evaluation is not possible, but a tabular overview is possible and can be used by the analyst to assess which solutions is the best to him or her.

Maryam Rahmani says:

Hello Mr. Linnart Felkl,
Thank you for sharing this approach. I
try to solve a problem in pulp.
I have a matrix in my problem. the dimension of this matrix is N*N*F*T.
That means the solver must calculate the binary N*N matrix for F times in T timeslot.
I use this matrix in my code with this method:

Afij_t = pulp.LpVariable.dict(”Afij_t,”,((t,f,i,j) for t in range(T) for f in range(F)
for i in range(N)for j in range(N)), cat =”binary”).

now i want result of this matrix in each timeslot becouse in each timeslot value of some of parameters and variables changing. with this method solver give me sum of all
time slot while i dont want it. i want cost of objective function in each time slot.
in other way i used N*N matrix for F times in a list with timeslot length but pulp calculate objective function costs for one timeslot and other timeslot ignore by pulp.
i want to know is ture this method for defintion matrix in pulp? if no what method
is true. sorry if my english language is not good. i hop that understand my intent.
thanks in advance.

Dear Sai,

I think this is already a very specific question. We would be happy to answer that as part of a paid consultation. For this you can write to me under below link, by filling out the contact form. I will direct you to the correct expert, who can answer your questions in future and propose a rate (i.e. some hourly consulting fee) to you.

Here is the contact form:
https://www.supplychaindataanalytics.com/contact/

RAHEEL SHAIKH says:

Hi, Thank you for the article.
I am new to optimization and PULP, can you guide if the optimization can be done where instead of constant1 x variable1 + constant2 x variable2 (ax + by), we have series of numbers i.e. a & b are series.

Cheers

Maryam Rahmani says:

Hello Linnart,

Thanks for share this approach in pulp.
i am trying to solve my problem with pulp. i have binary matrix like :

(Aij(f))(t)

that i and j have range(len(N)).
that means dimension of matrix is:
N*N*F*T.
a binary matrix N*N for F number execute in T timeslot.

i try to use this matrix in my code with this method:

Afij_t = pulp. LpVariable. dict(”Afij_t, ”, ((t, f, i, j) for t in range(T) for f in range(F) for i in range(N) for j in range(N)), cat =”binary”)

now i want to have result of solve
of this problem in each timeslots. but the solver gives results for sum of all timeslots. becouse my objective function must execute sum of all F number matrix (N*N) for T timeslot to give me best result in each timeslot becouse in each timeslot some constraints changing. also i try to use matrix with three dimension in a list with timeslots length but objective function cost execute just for one timeslot and another index of list ignore to solve by pulp.
i want to know is any way to access to each timeslots. sorry if my english languge is not good. i hop that you understand my intent. and so thanks for your comments and answers that they are so helpful.

Dear Maryam,

Thank you very much for your question. As I also wrote in another comment, I would like to forward such specific questions through our contact form, and have such topics answered as part of a paid consultation. The fee can be low, doesn`t need to be expensive. But as it is a lot of effort to support this blog I would like to forward this question to one of our experts. He or she will then get in touch, provide a rate (consulting fee estimate, e.g. on hourly basis), and will be able to answer any of your questions in future in this regard.

For this, please contact us via the contact form linked below:
https://www.supplychaindataanalytics.com/contact/

Monish Kumar says:

Hi Linnart,
Can you help me how to get the values of 2 objectives seperately in weighted method.

After solving the weighted problem, i.e. the problem with the scalarized objective function, you use its optimal solution to recalculate objective 1 and objective 2.

akashay says:

Hi,
how to take the print of the iterations? unbale to get using solutionTable.head()

Why is solutionTable.head() not working? Does solutionTable contain data or are you receiving an empty reference?

Shaifali says:

Hi,
How do i assume weights if there are more than 2 objectives?

For example by a) normalizing the objective function and b) applying a “3D” mesh of weight values. Its not difficult. Send me an e-mail and send me your problem and I can do it for you in 15 or 30min

Shaifali says:

Hi
Thank you so much for replying, so I am using a weighted sum objective which has 3 components say Obj1, Obj2, Obj3, right now I am assigning the weights manually, but now I want the code to decide optimal weights for all three objectives, I saw your example where you are using alpha and 1-alpha and giving a step size of 0.1 , i want to perform same on 3 objective problem which i have, i am using pulp as well

Shaifali says:

Hi
Thank you so much for replying
Sorry i am bugging you again, but i cannot share script, can you please provide me a link or example so that i can apply it 🙁

Ignacio says:

Hi, Linnart

Thanks for the post. I am looking to solve a location facility problem that minimizes costs of transport and costs of setup on the one hand and co2 emissions on the other. I have a hard time finding a simple example of such a scenario. Do you recommend some literature or source?

can be modelled as a mixed integer problem and solved with coin-or. You can model it with PuLP. There are some integer and MIP examples on the blog. to find them go ahead and use the search function

Mansour Eddaly says:

I address to the multiobjective permutation flowshop scheduling problem using the multiobjective mixed integer linear programming. I need a method to find the optimal set of non dominated solutions. I’am trying to use python with cplex. Please, how could you help me?
Many thanks

Dear Mr Mansour Eddaly, Thank you for your question. I can help you by implementing the solution for you. Please send me your file with the code in it, or your problem description, and I can give you a quote for what I would charge your for delivering a solution to you. Kind regards

btw, in my solution I will use an open source solver, as I do not currently hold a commercial solver. If you want to use CPLEX you would have to change the solver to CPLEX in the solution provided by me in the form of a Python source code.

Leave a Reply

Your email address will not be published. Required fields are marked *

This site uses Akismet to reduce spam. Learn how your comment data is processed.