The Australian National University Semester 2, 2023
School of Computing Assignment 1
COMP4691/8691: Optimisation
Release Date. 31 July 2023
Due Date. 18 August 2023, 18:00 AEST1
Maximum credit. 100 marks
Submission: Written answers go in an answers.pdf file with a limit of 2 pages (excluding an optional
title page). Code implementation go in the respective files specified in the questions below. Submission
and grading will be done on GitLab.
Requirements: This assignment requires usage of python3 and the numpy, scipy and PuLP packages.
Collaboration: You are allowed to discuss ideas with your fellow peers. However, you must write your
own written solutions and code by yourself with no external help. This means you are not allowed to look
at the solutions of your peers’ assignments.
Other notes:
• Some basic knowledge of numpy is required for the implementation questions. You can refer to
https://numpy.org/doc/stable/user/quickstart.html for almost everything you need.
• In written answers, marks may be lost if spelling, grammatical or markdown formatting errors, or
other forms of poor expression, make the meaning of your answer unclear or ambiguous.
• In written answers, marks may also be lost for supplying information that is irrelevant to answering
the question.
• Note that each subquestion in this assignment requires a written component.
Note: For the coding questions, you may use the packages that are already imported into your
solutions but do not import any other package. There will be a significant penalty if you import
extra packages. This also includes if you imported them but did not use them.
1Also see https://comp.anu.edu.au/courses/comp4691/deliverables/
1
Section 1: Simplex (50 Total Points)
In this section you will implement the simplex algorithm. As you have explored in lectures, the sim-
plex algorithm is a method for solving linear programs (LPs). The original simplex algorithm involves
maintaining a tableau for finding entering and exiting variables to emulate moving between basic feasible
solutions (BFS) at each step of the algorithm. We will instead implement the revised simplex method
which is equivalent to the original algorithm but instead uses a basis of the constraint matrix. We refer
to the revised simplex method whenever we mention the simplex algorithm.
The simplex algorithm
Assume m < n. We will be working with linear programs of the form
min
x∈Rn
c⊤x
s.t. Ax = b
x ≥ 0
b ≥ 0
where A ∈ Rm×n, b ∈ Rm, c ∈ Rn. We recall that we can represent vertices of the feasible region with
a basis B consisting of size m corresponding to indices of the variables x. In each step of the simplex
algorithm, we partition A, c, and x into two parts corresponding to basic and non-basic variables:
A¯ ∈ Rm×m A˜ ∈ Rm×(n−m)
c¯ ∈ Rm c˜ ∈ Rn−m
x¯ ∈ Rm x˜ ∈ Rn−m.
BFS. We are then able to find the basic variables of the BFS corresponding to the basis B by solving
A¯x¯ = b. (1.1)
Then our BFS is given by x˜ = 0 and x¯ as the solution to this equation.
Entering variable. It is possible to derive that the rate of change λ of the objective relative to the
equality constraint for the variables corresponding to the basis can be computed by solving
A¯⊤λ = c¯. (1.2)
We then have that
µ = c˜− A˜⊤λ (1.3)
represents the relative increase in objective relative to each non-basic variable. Thus negative entries in
µ indicate that the corresponding non-basic variables will decrease the objective if it enters the basis. If
there are no negative entries, then we have achieved the optimal solution. Otherwise, we use a pivot rule
to choose the entering variable from variables corresponding to the negative entries in µ.
Exiting variable. Given an entering variable k, denote A:,k the k-th column of A. Then we have that
d, the decrease in each basic variable for each unit increase in the entering variable, is computed by
A¯d = A:,k. (1.4)
Thus,
ν =
x¯
d
(1.5)
is a vector representing the ‘time’ it takes for each basic variable to reach zero as we increase the value
of the entering variable. If no elements are both strictly positive and finite, then we have an unbounded
problem. Otherwise, we use a pivot rule to choose the exiting variable from variables corresponding to
variables in ν that are both strictly positive and finite.
2
Solving an LP
To solve a linear program, we actually require running two iterations of the simplex algorithm described
previously: once on a transformed problem with a known basis for finding a BFS for the original problem,
and once on the main problem by using the BFS found in the modified problem as a starting point. To
construct the transformed problem, we introduce slack variables for each equality constraint in order to
get a trivial BFS with the slack variables as basis (you have to find out what values they take):
min
x∈Rn,s∈Rm
1⊤s
s.t. Ax+ s = b
x ≥ 0
s ≥ 0
b ≥ 0
Once the problem is solved, if any slack variable is non-zero, this means that the original problem is
infeasible as there is no solution where the equality constraints hold. Otherwise, all the slack variables
will no longer be in the basis and thus we have a BFS for the original problem.
Question 1.1: Problem modification for finding a starting BFS (10 Points)
• Why do we assume b ≥ 0 in the equality plus non-negative form for the simplex algorithm
and is it possible to do so for all linear programs? Justify your answer.
• Why do we assumem < n in the same form and is it possible to do so for all linear programs?
Justify your answer.
• Complete the implementation of the problem modification for finding a starting BFS in the
starting bfs lp function in the simplex.py file.
Hint: you will not be able to test your implementation for this subquestion yet until you finish
the next subquestion.
3
Question 1.2: The simplex algorithm (30 Points)
• Implement the simplex algorithm using the provided framework in the simplex and
simplex step functions with a pivot rule of your choice for the entering and exiting variables.
• Report and describe the pivot rules that you have chosen to implement and make a comment
on their advantages and disadvantages.
You should put your implementation of the simplex algorithm iterations into the simplex step
function. You may use additional helper functions if you feel necessary. You must not modify the
input and output types of these functions. To run your code, execute the command
python3 code/simplex.py data/q1-lp-x.json
for any x ∈ {a, b, c}.
Notes and hints:
• You should also handle the case where the problem is unbounded or infeasible even though
these are not provided in the test problems.
• Marks are awarded for the correctness and efficiency of the implementation.
• Code must be well formatted and commented.
• It may be a good idea to test your simplex algorithm implementation against a PuLP solver.
Optimising the simplex algorithm
Each step of the simplex algorithm requires solving a system of equations three times, as in (1.1), (1.2)
and (1.4), which is quite costly. However, an optimisation we can perform is to notice that the coefficient
matrix is the same in all three equations which means we can use LU factorisation to speed up the solving
process.
Question 1.3: LU factorisation optimisation (10 Points)
• Implement the LU factorisation optimisation in the simplex step opt function.
• Run experiments with both simplex algorithm implementations and report your results.
Explain all of your observations, whether they were expected or unexpected.
To run your code with the LU factorisation optimisation, add the --opt flag to the previous
command.
Hint: you can refer to https://en.wikipedia.org/wiki/LU_decomposition to understand LU
factorisation and how you may use it to solve linear equations. The lu factor and lu solve
functions from the scipy.linalg package may be very useful.
4
Section 2: Modelling (50 Total Points)
In this section you will model and clear an electricity market, where generators bid to supply power, with
a linear program. We will give you the description of the problem from which you will derive your own
variables, linear objective and constraints for modelling the problem.
The Market Problem
The market spans T consecutive time steps 0, . . . , T−1 and will be jointly cleared for a set of G generators.
Each time step has a fixed amount of demand that must be met by supply. Each generator has submitted
zero, one or more bids to supply power for each time step, where each bid consists of a price and
quantity. For example, a generator bid of (50, 100) indicates that the generator is willing to supply
up to 100MW where they get paid at least 50$ per MW for doing so. We also note that the market
may partially accept a bid. Continuing the previous example, the market is able to select a supply value
anywhere between 0 and 100MW for that bid. The cost associated with a (partial) bid is equal to the bid
price times the accepted supply for that bid.
We will also consider the generator ramp limit constraints. The limits constrain how much a genera-
tor can change their overall supply between two consecutive time steps indicated by max ramp up and
max ramp down for each generator. These values are always nonnegative. If a generator has multiple bids
for a particular time step, the overall supply for that generator in that time step is the sum of its accepted
supply across those bids. For the sake of this problem, generators can supply any amount of power in the
0-th time step.
The objective of the market is to accept the minimise the sum of bid costs under the above two sets of
constraints (supply equals demand for each time step, and generator ramp limits are followed).
Question 2.1: Model and solve the market problem (35 Points)
• Use an LP to implement the clear market function in the modelling.py file to model and
solve a MarketProblem problem instance described in The Market Problem section above by
returning a MarketSolution. You should now use a PuLP solver instead of your own solver.
• Report (1) the objective value of the solved linear program and (2) the highest bid price
that was accepted, either fully or partially, and the associated time steps and generators.
To run your code, execute the command
python3 code/modelling.py data/q2-bids-energy.json
To test your implementation on a smaller instance, add the --test flag to the previous command.
Notes and hints:
• Marks are awarded for the correctness and efficiency of the implementation.
• Your code should be able to run within 60 seconds. lpSum and lpDot may be helpful. As a
guideline, the tutors’ solutions run in a few seconds.
• Code must be well formatted and commented.
• Formulate the problem mathematically by hand before coding it up. You do not need to
report your mathematical formulation.
• Some bid prices are negative which represents the willingness of generators to pay the market
to supply a certain quantity. This may seem counter-intuitive but it is sometimes more cost
effective to pay to sell energy and keep a generator running than to turn a generator off.
• A real life application would use units in terms of energy (MWh) rather than power (MW).
We use MW to simplify the problem.
5
After clearing the market, generators are paid for their accepted supply at the marginal price for each
time step. The marginal price is a dual variable of your linear program associated with the constraint
that enforces the supply is equal to the demand in each time step.
Question 2.2: Dual variables (15 Points)
• Provide the marginal prices in the marginal prices field of the MarketSolution solution.
• Produce and report a timeseries plot of the marginal prices. This requires you to write your
own code in visualise solution.
• Explain how the marginal prices relate to the accepted bid prices and whether generators
would prefer being paid at their accepted bid prices or these marginal prices.
• Hypothesise and explain what would be the relationship between accepted bid prices and
marginal prices if we did not implement the generator ramp limit constraints.
Hint: Dual variables depend on how constraints are written and depend on the convention used
by the solver. Thus, you may get a result that is the negative of your intended quantity.