goalp is an R package for linear goal programming. It
allows solving basic, weighted, and lexicographic linear goal
programming problems, as well as a mixture of the weighted and
lexicographic approaches. The package accepts a human-readable
representation of the problem, and it automatically sets up the
corresponding linear programming problem, which it later solves using
the lp_solve linear
solver, through its R interface lpSolve. The
implementation of goal programming is based on the notation and approach
described in Ignizio (1983).
goalpImagine you want to grow some vegetables and herbs in your garden. You can grow aubergine (A), broccoli (B), carrots (C), dill (D) or endive (E). Each vegetable requires a certain amount of land to grow per seedling you plant, as well as a certain amount of water and fertilizer. After a bit of research, you find out that each seedling of each vegetable needs the following amount of resources to grow.
| Aubergine | Broccoli | Carrot | Dill | Endive | |
|---|---|---|---|---|---|
| Land | 60 | 30 | 20 | 10 | 20 |
| Water | 0.8 | 0.2 | 0.3 | 0.1 | 0.3 |
| Fertiliser | 30 | 20 | 10 | 5 | 20 |
Where land is expressed in cm of a planting line, water in litres per week, and fertiliser in grams per season. These values are only meant as an example, and do not reflect real needs of these plants.
Now imagine you have prepared the land already, and you have 3,000 cm of planting line available, you have 50 litres of water available during a week, and you have bought 2,000 grams of fertiliser.
You would like to decide how much to plant of each vegetable in such a way that your land, water, and fertiliser consumption is as close as possible to your available resources. You can go a bit above or below each amount, as you can always extend the planting lines, get additional water from other sources, or buy extra fertiliser, or leave any amount of them unused. But you would like to minimise additional work or expenditure, as well as avoid any waste. How can you decide, then, what to plant while consuming as close as possible to the available resources?
Goal programming can help us when we have multiple restrictions (land, water and fertiliser in our example) from which we want to deviate as little as possible. To formulate our vegetable patch problem as a goal programming problem, we must first write our constraints, which we will call goals, as the following set of equations.
\[60A + 30B + 20C + 10D + 20E = 3000\] \[0.8A + 0.2B + 0.3C + 0.1D + 0.3E = 50\] \[30A + 20B + 10C + 5D + 20E = 2000\]
Additional to the goals listed above, we assume all decision variables (A, B, C, D, E) to be non-negative. This is an implicit constraints that will be included in all goal programming problems in goalp.
We are unlikely to find a solution that satisfies all goals perfectly, so we allow for negative and positive deviations in each goal. For a given possible solution (i.e. a set of values for the decision variable), the negative deviation measures by how much the left-hand side of the goal is lacking to reach the right-hand side. The positive deviation, on the other hand, measures by how much the left-hand side of the equation is above or in excess of the right-hand side. The deviations are always expressed as non-negative values.
For example, if we were to plant two dozens of each vegetable, i.e. \(A=B=C=D=E=24\), then the left hand side of the “land” goal would have a value of \(60*24 + 30*24 + 20*24 + 10*24 +20*24 = 3360\), which means our positive deviation would be equal to \(d_{land+}=3360-3000=360\), because we have an excess of 360 units on our left hand side, above our goal. On the other hand, if we evaluate the left-hand side of our water goal, we find it has a value of \(0.8*24 + 0.2*24 + 0.3*24 + 0.1*24 + 0.3*24 = 40.8\), meaning we are lacking 9.2 units to reach our goal, in other words \(d_{water-}=50 - 40.8 = 9.2\).
Finding values for A, B, C, D, and E that keep us as close as possible to our goals is equivalent to minimising the sum of all deviations. Therefore, we just need to solve the following linear programming problem:
\[Min_{A, B, C, D, E, d_{i-}, d_{i+}}
\sum_{i\in \{ land, water, fertiliser\} } d_{i-} + d_{i+}\]
Subject to the following set of constraints: \[60A + 30B + 20C + 10D + 20E + d_{land-} -
d_{land+}= 3000\] \[0.8A + 0.2B + 0.3C
+ 0.1D + 0.3E + d_{water-} - d_{water+} = 50\] \[30A + 20B + 10C + 5D + 20E + d_{fertiliser-} -
d_{fertiliser+} = 2000\] We can solve this problem using
goalp. To do it, we simply need to follow four simply
steps.
goalp package if it is not loaded
already.goalp function with our goals as an
argument.summary on our results to see them on the
screen.The process is exemplified in the following code snippet.
# Load library
library(goalp)
# Write goals to a text variable
goals <- "Land : 60*A + 30*B + 20*C + 10*D + 20*E = 3000
Water: .8*A + .2*B + .3*C + .1*D + .3*E = 50
Fertiliser: 30*A + 20*B + 10*C + 5*D + 20*E = 2000"
# Solve problem
gp <- goalp(goals)
#> Solving underlying linear optimisation problem.
# Print results to screen
summary(gp)
#> Goal programming model run by tradpa using goalp 0.3.1 on R 4.2.1 for
#> Windows
#>
#> Problem formulation:
#> Land : 60*A + 30*B + 20*C + 10*D + 20*E = 3000 | 1 1 | 1# 1#
#> Water : 0.8*A + 0.2*B + 0.3*C + 0.1*D + 0.3*E = 50 | 1 1 | 1# 1#
#> Fertiliser : 30*A + 20*B + 10*C + 5*D + 20*E = 2000 | 1 1 | 1# 1#
#>
#> Objective function value: 5
#>
#> Solution:
#> value type
#> A 0 integer
#> B 0 integer
#> C 100 integer
#> D 0 integer
#> E 50 integer
#>
#> Deviations:
#> d- d+ w- w+ p- p+
#> Land 0 0 1 1 1 1
#> Water 5 0 1 1 1 1
#> Fertiliser 0 0 1 1 1 1We observe that the solution is \(A=0, B=0, C=100, D=0, E=50\), i.e. planting 100 carrots and 50 endives. The summary also reports the value of the deviations, with \(d_{water-}=5\), meaning our solution wastes 5 litres of water (i.e. it only uses 45 litres).
A few things to note when using goalp:
goalp."2*x + 3*y = 10" is a valid goal, but
"2x + 3y = 10" is not.Now imagine the fertiliser has a very long shelf life, so we do not care about using less of it, as we can always use whatever is left during the next season. In formal terms, this means that we are not interested in the minimising \(d_{fertiliser-}\), or alternatively, that our fertiliser goal can be rewritten as follows.
\[30A + 20B + 10C + 5D + 20E <= 2000\] At the same time, imagine we really enjoy working in our plot of land, so we do not care about having to prepare more land for planting. In formal terms, this means that we are not interested in minimising \(d_{land+}\), or alternatively, that our land goal can be rewritten as follows.
\[60A + 30B + 20C + 10D + 20E >= 3000\]
We can solve the new problem using goalp using the
following code.
# Load library
library(goalp)
# Write goals to a text variable
goals <- "Land : 60*A + 30*B + 20*C + 10*D + 20*E >= 3000
Water: .8*A + .2*B + .3*C + .1*D + .3*E = 50
Fertiliser: 30*A + 20*B + 10*C + 5*D + 20*E <= 2000"
# Solve problem
gp <- goalp(goals)
#> Solving underlying linear optimisation problem.
# Print results to screen
summary(gp)
#> Goal programming model run by tradpa using goalp 0.3.1 on R 4.2.1 for
#> Windows
#>
#> Problem formulation:
#> Land : 60*A + 30*B + 20*C + 10*D + 20*E >= 3000 | 1 0 | 1# Inf#
#> Water : 0.8*A + 0.2*B + 0.3*C + 0.1*D + 0.3*E = 50 | 1 1 | 1# 1#
#> Fertiliser : 30*A + 20*B + 10*C + 5*D + 20*E <= 2000 | 0 1 | Inf# 1#
#>
#> Objective function value: 0
#>
#> Solution:
#> value type
#> A 0 integer
#> B 25 integer
#> C 150 integer
#> D 0 integer
#> E 0 integer
#>
#> Deviations:
#> d- d+ w- w+ p- p+
#> Land 0 750 1 0 1 Inf
#> Water 0 0 1 1 1 1
#> Fertiliser 0 0 0 1 Inf 1The solution now changes, with the optimal now being planting 25 broccoli and 150 carrots. This new solution requires preparing 750 additional cm of land (7.5 metres) for planting, but consumes all available water and fertiliser.
Now imagine we enjoy variety in our diet, so we would like to to produce at least 10 aubergines and 10 broccoli. While we like carrots, we do not want too many of them, so we would like to plant only between 10 and 20. Dill, on the other hand, we definitely do not need more than 5. We can add these ranges to the formulation of our problem.
To add lower, upper and both kinds of bounds to the decision
variables in our problem, we use the key words lBound,
uBound, and range respectively. Note that
these are hard constraints, as no deviations are allowed for them.
# Load library
library(goalp)
# Write goals to a text variable
goals <- "Land : 60*A + 30*B + 20*C + 10*D + 20*E >= 3000
Water: .8*A + .2*B + .3*C + .1*D + .3*E = 50
Fertiliser: 30*A + 20*B + 10*C + 5*D + 20*E <= 2000
A lBound 10
B lBound 10
C range [10,20]
D uBound 5"
# Solve problem
gp <- goalp(goals)
#> Solving underlying linear optimisation problem.
# Print results to screen
summary(gp)
#> Goal programming model run by tradpa using goalp 0.3.1 on R 4.2.1 for
#> Windows
#>
#> Problem formulation:
#> Land : 60*A + 30*B + 20*C + 10*D + 20*E >= 3000 | 1 0 | 1# Inf#
#> Water : 0.8*A + 0.2*B + 0.3*C + 0.1*D + 0.3*E = 50 | 1 1 | 1# 1#
#> Fertiliser : 30*A + 20*B + 10*C + 5*D + 20*E <= 2000 | 0 1 | Inf# 1#
#> A lBound 10
#> B lBound 10
#> C lBound 10
#> C uBound 20
#> D uBound 5
#>
#> Objective function value: 0
#>
#> Solution:
#> value type
#> A 52 integer
#> B 12 integer
#> C 20 integer
#> D 0 integer
#> E 0 integer
#>
#> Deviations:
#> d- d+ w- w+ p- p+
#> Land 0 880 1 0 1 Inf
#> Water 0 0 1 1 1 1
#> Fertiliser 0 0 0 1 Inf 1
#> A_lBound NA 42 NA 0 NA Inf
#> B_lBound NA 2 NA 0 NA Inf
#> C_lBound NA 10 NA 0 NA Inf
#> C_uBound 0 NA 0 NA Inf NA
#> D_uBound 5 NA 0 NA Inf NAThe optimal solution is now planting 52 aubergines, 12 broccoli and 20 carrots. Note that this solution satisfy all bounds imposed on the problem. It also implies preparing 8.80 additional metres of land for planting.
Lower boundaries have the negative deviation \(d_{i-}\) set to NA, while
upper boundaries have the positive deviation \(d_{i+}\) set to NA. This is
because lower boundaries do not allow negative deviations, as the
variable must be equal or bigger than the lower boundary. Similarly,
upper bounds do not allow for positive deviations, therefore their
values are shown as NA.
Note that, unlike goals, names cannot be assigned to bounds.
So far, we have given the same weight to all deviations. For example, lacking 1 cm of land was just as bad as using one litre more or less than 50 litres. Very often, this is not the case. Continuing with the example, deviating by one litre of water is probably more important than deviating by 1 cm of land.
To address this disparity, we can weight the deviations in the
objective function differently for each goal. For example, we could say
that what matters are the deviations in terms of percentage, so
deviating by 30 cm (1% of 3000 cm) is just as important as deviating by
0.5 litres (1% of 50 litres), or 20 grams of fertiliser (1% of 2000
grams). To do this, we need to assign weights to each deviation. We can
do this easily in goalp, as follows.
# Load library
library(goalp)
# Write goals to a text variable
goals <- "Land : 60*A + 30*B + 20*C + 10*D + 20*E >= 3000 | 1/3000
Water: .8*A + .2*B + .3*C + .1*D + .3*E = 50 | 1/50
Fertiliser: 30*A + 20*B + 10*C + 5*D + 20*E <= 2000 | 1/2000
A lBound 10
B lBound 10
C range [10,20]
D uBound 5"
# Solve problem
gp <- goalp(goals)
#> Solving underlying linear optimisation problem.
# Print results to screen
summary(gp)
#> Goal programming model run by tradpa using goalp 0.3.1 on R 4.2.1 for
#> Windows
#>
#> Problem formulation:
#> Land : 60*A + 30*B + 20*C + 10*D + 20*E >= 3000 | 0.0003333 0 | 1# Inf#
#> Water : 0.8*A + 0.2*B + 0.3*C + 0.1*D + 0.3*E = 50 | 0.02 0.02 | 1# 1#
#> Fertiliser : 30*A + 20*B + 10*C + 5*D + 20*E <= 2000 | 0 0.0005 | Inf#
#> 1#
#> A lBound 10
#> B lBound 10
#> C lBound 10
#> C uBound 20
#> D uBound 5
#>
#> Objective function value: 0
#>
#> Solution:
#> value type
#> A 52 integer
#> B 10 integer
#> C 20 integer
#> D 4 integer
#> E 0 integer
#>
#> Deviations:
#> d- d+ w- w+ p- p+
#> Land 0 860 0.0003333333 0e+00 1 Inf
#> Water 0 0 0.0200000000 2e-02 1 1
#> Fertiliser 20 0 0.0000000000 5e-04 Inf 1
#> A_lBound NA 42 NA 0e+00 NA Inf
#> B_lBound NA 0 NA 0e+00 NA Inf
#> C_lBound NA 10 NA 0e+00 NA Inf
#> C_uBound 0 NA 0.0000000000 NA Inf NA
#> D_uBound 1 NA 0.0000000000 NA Inf NAChanging the weights of each deviation leads to a different optimal solution, this time the solution being 52 aubergines, 10 broccoli, 20 carrots, and 4 dill.
The printed solution shows the weight of each deviation next to the
value of each deviation, under the columns w- and
w+ for the weights of the negative and positive deviations,
respectively.
Note that goals expressed as equalities (=) have weights
different from zero for both the negative (\(d_{i-}\)) and positive deviations (\(d_{i+}\)). Goals expressed as “bigger than”
inequalities (>=) have the weight of \(d_{i+}\) automatically set to zero, as
going over the goal is not penalised. Goals expressed as “smaller than”
inequalities (<=) have the weight of \(d_{i-}\) automatically set to zero, as
being under the goal is not penalised.
The same effect could be achieved by setting the normW
options to TRUE when calling goalp,
i.e. goalp(goal, normW=TRUE). The option normW
stands for “normalised weights”. It automatically defines weights as the
inverse of the right-hand side of each goal.
Note that bounds do not accept the use of weights.
Now imagine that not only we want to reduce deviations from our ideal water usage of 50 litres, but also it is more onerous for us to exceed 50 litres, than it is to use less than that. For example, it could be that additional water has to be fetched from a distant tap, so it requires a significant effort, while using less water is of little annoyance.
We can take this into consideration by defining different weights for the negative and positive deviations for the “water” goal. We do this by defining two weights for the goal, first the negative deviation weight, followed by the weight of the positive deviation, as shown in the example below.
# Load library
library(goalp)
# Write goals to a text variable
goals <- "Land : 60*A + 30*B + 20*C + 10*D + 20*E >= 3000 | 1/3000
Water: .8*A + .2*B + .3*C + .1*D + .3*E = 50 | 1/50 10/50
Fertiliser: 30*A + 20*B + 10*C + 5*D + 20*E <= 2000 | 1/2000
A lBound 10
B lBound 10
C range [10,20]
D uBound 5"
# Solve problem
gp <- goalp(goals)
#> Solving underlying linear optimisation problem.
# Print results to screen
summary(gp)
#> Goal programming model run by tradpa using goalp 0.3.1 on R 4.2.1 for
#> Windows
#>
#> Problem formulation:
#> Land : 60*A + 30*B + 20*C + 10*D + 20*E >= 3000 | 0.0003333 0 | 1# Inf#
#> Water : 0.8*A + 0.2*B + 0.3*C + 0.1*D + 0.3*E = 50 | 0.02 0.2 | 1# 1#
#> Fertiliser : 30*A + 20*B + 10*C + 5*D + 20*E <= 2000 | 0 0.0005 | Inf#
#> 1#
#> A lBound 10
#> B lBound 10
#> C lBound 10
#> C uBound 20
#> D uBound 5
#>
#> Objective function value: 0
#>
#> Solution:
#> value type
#> A 52 integer
#> B 10 integer
#> C 20 integer
#> D 4 integer
#> E 0 integer
#>
#> Deviations:
#> d- d+ w- w+ p- p+
#> Land 0 860 0.0003333333 0e+00 1 Inf
#> Water 0 0 0.0200000000 2e-01 1 1
#> Fertiliser 20 0 0.0000000000 5e-04 Inf 1
#> A_lBound NA 42 NA 0e+00 NA Inf
#> B_lBound NA 0 NA 0e+00 NA Inf
#> C_lBound NA 10 NA 0e+00 NA Inf
#> C_uBound 0 NA 0.0000000000 NA Inf NA
#> D_uBound 1 NA 0.0000000000 NA Inf NAIn this case the result does not change, as the previous solution had
negative and positive deviations equal to zero for the water goal.
Nevertheless, we can see that the weights reported in the
w- and w+ columns are different for the water
goal.
A more concise way to define the previous problem would have been the
following. As we use the normW=TRUE option, the scaling of
weights by the inverse of the righ-hand side value of each goal is
automatic.
# Load library
library(goalp)
# Write goals to a text variable
goals <- "Land : 60*A + 30*B + 20*C + 10*D + 20*E >= 3000
Water: .8*A + .2*B + .3*C + .1*D + .3*E = 50 | 1 10
Fertiliser: 30*A + 20*B + 10*C + 5*D + 20*E <= 2000
A lBound 10
B lBound 10
C range [10,20]
D uBound 5"
# Solve problem
gp <- goalp(goals, normW=TRUE)
# Print results to screen
summary(gp)Now imagine that we are more interested in using all of our available land, than in fulfilling our water and fertiliser goals. For example, it could be that all the land we do not use is assigned to other people in the next harvest season. When one or more goals are more important than others, we use lexicographic goal programming.
Lexicographic goal programming requires solving the problem multiple times, but considering only some deviation in each occasion. The idea idea is that we first optimise for the most important goal (or deviation), we fix those deviations to the best we found, and then we optimise again considering only the next most important deviation, and so on and so forth.
For example, let’s say our first priority is using at least 3000 cm
of land, in other words, we want to minimise \(d_{land-}\), so we assign it a priority of
1. Our next priority is avoiding to use more than 50 litres of water, so
we assign \(d_{water+}\) a priority of
2. Then we care the same about all remaining deviations, so we assign
priority 3 to each of them. We communicate this priority to
goalp using the following syntax.
# Load library
library(goalp)
# Write goals to a text variable
goals <- "Land : 60*A + 30*B + 20*C + 10*D + 20*E >= 3000 | 1#
Water: .8*A + .2*B + .3*C + .1*D + .3*E = 50 | 3# 2#
Fertiliser: 30*A + 20*B + 10*C + 5*D + 20*E <= 2000 | 3#
A lBound 10
B lBound 10
C range [10,20]
D uBound 5"
# Solve problem
gp <- goalp(goals)
#> Solving for priority level 1 (1/3).
#> Solving for priority level 2 (2/3).
#> Solving for priority level 3 (3/3).
# Print results to screen
summary(gp)
#> Goal programming model run by tradpa using goalp 0.3.1 on R 4.2.1 for
#> Windows
#>
#> Problem formulation:
#> Land : 60*A + 30*B + 20*C + 10*D + 20*E >= 3000 | 1 0 | 1# Inf#
#> Water : 0.8*A + 0.2*B + 0.3*C + 0.1*D + 0.3*E = 50 | 1 1 | 3# 2#
#> Fertiliser : 30*A + 20*B + 10*C + 5*D + 20*E <= 2000 | 0 1 | Inf# 3#
#> A lBound 10
#> B lBound 10
#> C lBound 10
#> C uBound 20
#> D uBound 5
#>
#> Objective function value: 0
#>
#> Solution:
#> value type
#> A 52 integer
#> B 12 integer
#> C 20 integer
#> D 0 integer
#> E 0 integer
#>
#> Deviations:
#> d- d+ w- w+ p- p+
#> Land 0 880 1 0 1 Inf
#> Water 0 0 1 1 3 2
#> Fertiliser 0 0 0 1 Inf 3
#> A_lBound NA 42 NA 0 NA Inf
#> B_lBound NA 2 NA 0 NA Inf
#> C_lBound NA 10 NA 0 NA Inf
#> C_uBound 0 NA 0 NA Inf NA
#> D_uBound 5 NA 0 NA Inf NAThe solution has now changed slightly to 52 aubergines, 12 broccoli,
and 20 carrots. As the previous solution already did not use less than
2000 cm of land the solution was not expected to change drastically.
When we run goalp, we see that three sub-problems are
solved, one for each level of priority.
Note that, unlike weights, the syntax for priority includes a hashtag
(#).
It is also possible to combine lexicographic priorities and weights.
But then the weights will be relevant only within the same level of
priority. We could, for example, obtain sligtly different results if we
instead run goalp(goals, normW=TRUE).
When printing the results, the lexicographic priority of each
deviation is shown next to its weights, in columns p- and
p+ in the table at the bottom of the summary.
goalp by default assumes that all decision variables are
integers. But they can also be defined as continuous
(non-negative real numbers) or binary (variables that can only
take values 0 or 1). The typre of each variable can be specified through
the varType option when calling goalp, for
example
goalp(goals, varType=c(x="cont", y="int", z="bin")), where
x, y and z are the names of the
decision variables.
goalp allows for enforced equality constraints using the
== symbol. This will turn a goal into an effective
constraint. Note that this increases the chances of the solver not
finding a solution for the underlying linear optimisation problem.
goalpalso allows for enforcing of constraints by setting goal's weights equal toNA. If a weight is set toNA`,
then the corresponding deviation is removed from the underlying linear
optimisation problem. This is specially useful in the case of enforcing
inequalities. For example, to enforce \(x + y
>= 5\):
goals <- "g1: 3*x + 2*y + z = 20
g2: x + y >= 5 | NA"
gp <- goalp(goals)
summary(gp)Besides the text syntax describing a goal programming problem,
goalp also allows defining a problem using matrices. This
is specially useful in the case of large problems. The problem is
expressed in the form Ax = b, with the user having to
provide the coefficient matrix A, with as many rows as
goals and columns as decision variables, the vector b, with
the right-hand side values of the goals, and an additional character
vector m indicating if each goal relates the left- and
right-hand sides using "=", ">=",
"<=", "==", "lBound", or
"uBound". Type ?goalp in the console for more
details.
Ignizio, J.P (1983) Generalized goal programming: An overview. Computers and Operations Research 10, 277-289. doi