Report Submission Information (must be completed before submitting report!)
• Student 1 Full Name and Number : Hao Xu 718329
• Student 2 Full Name and Number :
• Workshop day : Thursday
• Workshop time : 6:15 – 8:15pm
Workshop 1 – Optimisation [2 weeks] ¶
Objectives:¶
• Learn how to formulate optimisation problems in practice.
• Familiarise yourself with practical software tools used for optimisation.
• Solve optimisation problems using Python Scipy and Matlab.
• Connect theoretical knowledge and practical usage by doing it yourself.
Common objectives of all workshops: Gain hands-on experience and learn by doing! Understand how theoretical knowledge discussed in lectures relates to practice. Develop motivation for gaining further theoretical and practical knowledge beyond the subject material.
Overview:¶
Another name for the field of “Optimisation” is “Mathematical Optimisation.” As the name indicates optimisation is an area of applied mathematics. It is possible to study optimisation entirely from a mathematical perspective. However, engineers are interested in solving real-world problems in a principled way. Many engineering problems can be and are formulated as optimisation problems. In those cases, mathematical optimisation provides a solid theoretical foundation for solving them in a principled way.
In this workshop, you will learn how to formulate and solve optimisation problems in practice. This will give you a chance to connect theoretical knowledge and practical usage by doing it yourself. You will familiarise yourself with practical optimisation tools for Python. These are chosen completely for educational reasons (simplicity, accessibility, cost). While the underlying mathematics is timeless, optimisation software evolves with time, and can be diverse. Fortunately, once you learn one or two, it should be rather easy to learn others now and in the future, because software designers often try to make it user friendly and take into account what people already know.
In the future, you should consider and learn serious optimisation software for scalability and reliability. They can be complex and/or expensive but they get the job done for serious engineering. Teaching such software takes too much time and is beyond the scope of this subject.
Workshop Preparation: [before you arrive to the lab]¶
You can come to the workshops as you are or you can prepare beforehand to learn much more! We will give you a lot of time to finish the tasks but those are the bare minimums. Just like in the lectures, the topics we cover in the workshops are quite deep and we can only do so much in two hours. There is much more to learn and coming prepared to the workshop is one of the best ways to gain more knowledge! For example, there are a few questions in each workshop which you can answer beforehand.
Self-learning is one of the most important skills that you should acquire as a student. Today, self-learning is much easier than it used to be thanks to a plethora of online resources. For this workshop, start by exploring the resource mentioned in the preparation steps below.
Workshop Preparation Steps:¶
1. Common step for all workshops: read the Workshop Manual (Jupyter Notebook) beforehand!
2. Review relevant lecture slides on optimisation.
3. Read/check relevant reading material and links from LMS/Resources-Reading
4. Check the embedded links below hints and background.
5. [optional] You can start with workshop tasks and questions
Tasks and Questions:¶
Follow the procedures described below, perform the given tasks and answer the workshop questions on the Python notebook itself! The marks associated with each question are clearly stated. Keep your answers to the point and complete to get full marks! Ensure that your code is clean and appropriately commented.
The resulting notebook will be your Workshop Report!
The goal is to learn, NOT blindly follow the procedures in the fastest possible way! Do not simply copy-paste answers (from Internet, friends, etc.). You can and should use all available resources but only to develop your own understanding. If you copy-paste, you will pay the price in the final exam!
Section 1: Convex Functions¶
Remember the definition of convex and concave functions from lecture slides. Functions are mathematical objects but they are used in engineering in very practical ways, for example, to represent the relationship between two quantities. Let’s draw a function!
In [1]:
%matplotlib notebook
import numpy as np
import matplotlib.pyplot as plt
# define function f(x)
def f(x):
return 3*(x-1)**2
# define x and y
x = np.linspace(-20, 20, 100) # 100 equally spaced points on interval [-20,20]
y = f(x) # call function f(x) and set y to the function’s return value
# Plot the function y=f(x)
plt.plot(x,y)
plt.xlabel(‘x’)
plt.ylabel(‘y’)
plt.title(‘$y=3(x-1)^2$’)
plt.show()

Question 1.1. (1 pt)¶
Plot one concave and one nonconvex functions of your choosing (preferably in 3D). Provide their formulas below.
In [50]:
import numpy as np
import matplotlib.pyplot as plt
#concave function
def f1(x):
return -(x-1)**2
#nonconvex function
def f2(x):
return (np.sin(x))
# define x and y
x = np.linspace(-20, 20, 100) # 100 equally spaced points on interval [-20,20]
y1 = f1(x) # call function f(x) and set y to the function’s return value
y2 = f2(x)
# Plot the function y=f(x)
plt.subplot(2,1,1)
plt.plot(x,y1)
plt.xlabel(‘x’)
plt.ylabel(‘y1’)
plt.title(‘$y1=-(x-1)^2$’)
plt.subplot(2,1,2)
plt.plot(x,y2)
plt.xlabel(‘x’)
plt.ylabel(‘y2’)
plt.title(‘$y2=sin(x)$’)
plt.subplots_adjust(hspace =1)
plt.show()

Answer as text here 1st Function: y1=-(x-1)^2 2nd Function: y2=sin(x)
Question 1.2. (1 pts)¶
How would you determine whether a single or multi-variate continuously differentiable function is convex or not?
Note that the question becomes very tricky if you have a parametric multivariate polynomials of degree four or higher!
[Optional] An interesting paper (for those who wish to go deeper) http://web.mit.edu/~a_a_a/Public/Publications/convexity_nphard.pdf
The function is convex if f(ax+(1-a)y)<= af(x) + (1-a)f(y) for x,y are the two points from the function set and a is in [0,1]
Question 1.3. (2 pts)¶
Why are convex optimisation problems considered to be easy to solve? Consider optimality conditions of unconstrained functions in your answer. Plot first and second order derivative functions of the concave and non-convex functions you have chosen above (as part of Question 1.1) to further support your argument.
Considering the concave function as y1, when finding the global maximum, it is similar to the convex. By looking at the first order derivative. It is clearly to see the global maximum is when x = 1 which means we can easily find global maximum in concave case, so as the global minimum in convex case, but for nonconvex question as y2, it has lots of zeros in the two derivative functions, and there is not clearly which zero point is the global maximum or minimum.
In [89]:
#Because the local minimum value of convex function is guaranteed to be the global minimum
#and convex problem are guaranteed to have globally optimal solutions
import numpy as np
import matplotlib.pyplot as plt
import pylab
#concave function
def f11(x):
return -2*x+2
#nonconvex function
def f12(x):
return -2*np.ones(len(x))
def f21(x):
return(np.cos(x))
def f22(x):
return(np.sin(-x))
# define x and y
x = np.linspace(-20, 20, 100) # 100 equally spaced points on interval [-20,20]
y11 = f11(x)
y12 = f12(x)
y21 = f21(x)
y22 = f22(x)
# Plot the function y=f(x)
plt.figure(figsize = (10,20))
plt.subplot(4,1,1)
plt.plot(x,y11)
plt.ylabel('y11')
plt.title('$y1\'=-2*x+2$')
plt.subplot(4,1,2)
plt.plot(x,y12)
plt.ylabel('y12')
plt.title('$y1\'\'=-2$')
plt.subplot(4,1,3)
plt.plot(x,y21)
plt.ylabel('y21')
plt.title('$y2\'=cos(x)$')
plt.subplot(4,1,4)
plt.plot(x,y22)
plt.xlabel('x')
plt.ylabel('y22')
plt.title('$y2\'\'=-sin(x)$')
plt.subplots_adjust(hspace =1)
plt.show()

Section 2: Unconstrained Optimisation¶
2.1 Example. Aloha communication protocol¶

Aloha is a well-known random access or MAC (Media/multiple Access Control) communication protocol. It enables multiple nodes sharing a broadcast channel without any additional signaling in a distributed manner. Unlike FDMA or TDMA (frequency or time-division multiple access), the channel is not divided into segments beforehand and collisions of packets due to simultaneous transmissions by nodes are allowed. In slotted Aloha, the nodes can only transmit at the beginning of time slots, which are kept by a global/shared clock. See Aloha for further background information.
Slotted Aloha Efficiency¶
For an $N$-node slotted Aloha system, where each node transmits with a probability $p$, the throughput of the system is given by $$ S(p) = N p (1-p)^{N-1}$$
Question 2.1. (2 pts)¶
Formally define the optimisation problem to find the optimal probability $p$ that maximises the throughput. Clearly identify the objective and decision variable(s). Is the objective convex or concave? Show through derivation.
Note that, there is the constraint $0 \leq p \leq 1$ on probability $p$ but we will ignore it for now.
First we conclude the first derivative of the S, $\frac{ds}{dp} = N * {(1-p)}^{(N-2)}*[1-p-p*(N-1)]$, so the zero point at $1-N*p = 0$, which lead to $p = \frac{1}{N}$. This is the optimal probability p that maximises the throughput S. By taking two points around 1/N which are $\frac{1}{2N}$ and $\frac{2}{N}$ respectively, the first derivative are $\frac{1}{2}$ and -1, which states that the objective is concave.
Question 2.2. (2 pts)¶
Plot the performance function and its derivative for $N=10$ nodes. Is the objective function convex or concave? Determine using mathematical methods. Investigate the property of the derivative function of the objective. What do you call such functions?
Answer as text here
In [2]:
import numpy as np
import matplotlib.pyplot as plt
import pylab
import math
N = 10;
def f(x):
return N*x*(1-x)**(N-1)
x = np.linspace(0,1,100)
y = f(x)
plt.plot(x,y)
plt.xlabel('p')
plt.ylabel('S(p)')
plt.title('${S(p) = N * p * (1-p)}^{(N-1)}$')
plt.show()

Question 2.3. (2 pts)¶
Find the optimal probability $p$ for $N=10$ nodes. Use an appropriate package from Scipy. Cross-check your answer with a mathematical formula that you should derive by hand.
In [39]:
from scipy import optimize as opt
import numpy as np
import matplotlib.pyplot as plt
import pylab
import math
N = 10;
def f(x):
return -N*x*(1-x)**(N-1)
x = np.linspace(0,1,100)
x0 = 0
y = f(x)
pmax = opt.minimize(f,x0)
print('At N = 10, Pmax =' ,pmax.x[0])
At N = 10, Pmax = 0.10000002145833721
As the function is concave and there is only minimize function in the scipy package, we need to convert the concave function to a convex one first and then minimize the convex function which will give us the result of the maximum point the concave function. As the previous statements, the result should be $p = \frac{1}{N}$ and $N = 10$ which is consistent with the result in Question2.3
2.2 Gradient Descent Algorithms¶
Question 2.4. (10 pts)¶
Write your own gradient algorithm (with constant step size) to solve the problem $$ \max_x x^T Q x + r^T x,$$ where $x,\, r \in \mathbb{R}^2$ and $Q$ is a $2\times 2$ positive definite matrix of your choice. Cross-check your answers using one of the standard optimisation packages, e.g. scipy or cvxpy.
1. If Q is positive definite, then what type of optimisation problem is this? Give your answer using mathematical tools learned in classroom. Would anything change if $Q$ were not positive definite? Plot both cases and comment.
2. Focusing on positive definite $Q$, what happens if you choose your fixed step size too large or too small? Observe and comment.
3. Choose $Q$ in such a way that is has a low [condition number] (https://www.encyclopediaofmath.org/index.php/Condition_number), see also. Next, choose a $Q$ with a high condition number. Compare the performance of your algorithm in both cases, and comment.
4. Now, solve both versions of the problem using (b) diminishing step size (c) Armijo rule/wolf test (line search). Discuss stopping criteria for all variants.
5. Plot your trajectories clearly showing the iterations of the gradient algorithm. You should also show either by displaying the level sets of the objective or the objective function itself (if resorting to a 3D plot).
Answer as text here Solution at: https://en.wikipedia.org/wiki/Gradient_descent
In [6]:
''' Answer as code here '''
Out[6]:
' Answer as code here '
Question 2.5. (2 pts)¶
Choose one of the versions of the problem and algorithms in Question 2.4 that finds the correct solution.
1. Use objective function itself or a norm of its gradient as a descent (Lyapunov) function that establishes the convergence of the solution algorithm $A(x)$. Plot the value of the chosen descent function versus the trajectory or steps.
2. Calculate and plot $||x(n)-x^*||$ over iterations $n=0,\ldots$, to establish that your solution algorithm leads to a pseudo-contraction.
Answer as text here
In [7]:
''' Answer as code here '''
Out[7]:
' Answer as code here '
Section 3: Constrained Optimisation¶
3.1 Example. Waterfilling in Communications¶
by Robert Gowers, Roger Hill, Sami Al-Izzi, Timothy Pollington and Keith Briggs. From the book by Boyd and Vandenberghe, Convex Optimization, Example 5.2 page 245.
$$\min_{x} \sum_{i=1}^N -\log(\alpha_i + x_i)$$
$$\text{subject to } x_i \geq 0, \; \forall i, \text{ and } \sum_{i=1}^N x_i = P $$
This problem arises in information/communication theory, in allocating power to a set of $n$ communication channels. The variable $x_i$ represents the transmitter power allocated to the i-th channel, and $\log(\alpha_i + x_i)$ gives the capacity or communication rate of the channel, where $\alpha_i>0$ represents the floor above the baseline at which power can be added to the channel. The problem is to allocate a total power of one to the channels, in order to maximize the total communication rate.
This can be solved using a classic water filling algorithm.

Question 3.1. (4 pts)¶
1. (1 pt) Is the problem in Example 3.1 convex? Formally explain/argue why or why not. What does this imply regarding the solution?
2. (2 pts) Solve the problem above for $N=8$ and a randomly chosen $\alpha$ vector. You can use for example Cvxpy package. Cross-check your answer with another software (package), e.g. Matlab or Scipy.
3. (1 pt) Write the Lagrangian, KKT conditions, and find numerically the Lagrange multipliers associated with the solution (using the software package/function). Which constraints are active? Explain and discuss briefly.
Answer as text here
Solution at https://www.cvxpy.org/examples/applications/water_filling_BVex5.2.html
In [8]:
”’ Answer as code here ”’
Out[8]:
‘ Answer as code here ‘
3.2 Example. Economic Dispatch in Power Generation¶
The problem is formulated as $$ \min_P \sum_{i=1}^N c_i P_i $$ $$\text{subject to } P_{i,max} \geq P_i \geq 0, \; \forall i, \text{ and } \sum_{i=1}^N P_i = P_{demand} $$
Here, $P_1,\ldots P_N$ are the power generated by Generators $1,\ldots,N$, $c_i$ is the per-unit generation cost of the i-th generator, and $P_{demand}$ is the instantaneous power demand that needs to be satisfied by aggregate generation. More complex formulations take into account transmission, generator ramp-up and down constraints, and reactive power among other things.
Question 3.2. (4 pts)¶
Let us get inspired from generation in Victoria with $N=12$ biggest generators that have more than 200MW capacity. Choose their maximum generation randomly or from the Victoria generator report if you wish to be more realistic. Generate a random cost vector $c$ varying between $10-50$ AUD per MWh. (Optionally, you can search and find how much different generation types cost if you are interested). Let the demand be $P_{demand}=5000MW$.
Solve this simplified economic dispatch problem defined above. The resulting merit order is the generation that would have been if there was no NEM (electricity market).
More about electricity market and generation at https://www.aemo.com.au/ See also this NEM overview introductory document (right click to download) and the Victoria generator report as of January 2019.
1. Solve the problem using cvxpy.
2. What type of an optimisation problem is this? Briefly explain.
3. Formulate by hand the dual problem and solve it with cvxpy. Is there a duality gap? Explain briefly why or why not.
4. Briefly comment on simplex algorithm and solve the problem using Scipy and simplex algorithm. Compare your results.
Answer as text here
In [9]:
”’ Answer as code here ”’
Out[9]:
‘ Answer as code here ‘
3.3 Example. Power Control in Wireless Communication¶
Adapted from Boyd, Kim, Vandenberghe, and Hassibi, “A Tutorial on Geometric Programming.”
The power control problem in wireless communications aims to minimise the total transmitter power available across $N$ trasmitters while concurrently achieving good (or a pre-defined minimum) performance.
The technical setup is as follows. Each transmitter $i$ transmits with a power level $P_i$ bounded below and above by a minimum and maximum level. The power of the signal received from transmitter $j$ at receiver $i$ is $G_{ij} P_{j}$, where $G_{ij} > 0$ represents the path gain (often loss) from transmitter $j$ to receiver $i$. The signal power at the intended receiver $i$ is $G_{ii} P_i$, and the interference power at receiver $i$ from other transmitters is given by $\sum_{k \neq i} G_{ik}P_k$. The (background) noise power at receiver $i$ is $\sigma_i$. Thus, the Signal to Interference and Noise Ratio (SINR) of the $i$th receiver-transmitter pair is
$$ S_i = \frac{G_{ii}P_i}{\sum_{k \neq i} G_{ik}P_k + \sigma_i }. $$
The minimum SINR represents a performance lower bound for this system, $S^{\text min}$.
The resulting optimisation problem is formulated as
$$ \begin{array}{ll} \min_{P} & \sum_{i=1}^N P_i \\ \text{subject to} & P^{min} \leq P_i \leq P^{max}, \; \forall i \\ & \dfrac{G_{ii}P_i}{\sigma_i + \sum_{k \neq i} G_{ik}P_k} \geq S^{min} , \; \forall i \\ \end{array} $$
Question 3.3. (10 pts)¶
Let $N=6$, $P^{min}=0.1$, $P^{max}=5$, $\sigma=0.2$ (same for all). Create a random path loss matrix $G$, where off-diagonal elements are between $0.1$ and $0.9$ and the diagonal elements are equal to $1$.
1. (2 pts) Write down the Langrangian and KKT conditions of this problem.
2. (2 pts) Solve the problem first with $S^{min}=0$ using cvxpy. Plot the power levels and SINRs that you obtain.
3. (2 pts) What happens if you choose an $S^{min}$ that is larger? Solve the problem again and document your results. What happens if you choose a very large $S^{min}$? Observe and comment.
4. (4 pts) Solve the problem using a combination of active set and penalty function methods. Specifically, choose a penalty function to impose the $S^{min}$ constraint and use an active set method for $P^{min}$, $P^{max}$ power constraints. Choose tighter constraints to make the problem more interesting.
Answer as text here
In [10]:
”’ Answer as code here ”’
Out[10]:
‘ Answer as code here ‘
3.4 (Optional Bonus, 10 pts) Model Predictive Control¶

It is possible to formulate a discrete-time, finite-horizon optimal-control as a constrained optimisation problem. This is quite useful since it allows making use of powerful optimisation solvers in addressing the control problem. This formulation is called Model Predictive Control (MPC).
Specifically, consider a system with a state vector $x_t\in {\bf R}^n$ that varies over the discrete time steps $t=0,\ldots,T$, and control actions $u_t\in {\bf R}^m$ that affect the state as part of a linear dynamical system formulated as
$$ x_{t+1} = A x_t + B u_t, $$
where $A \in {\bf R}^{n\times n}$ and $B \in {\bf R}^{n\times m}$ are system matrices.
The goal is to find the optimal actions $u_0,\ldots,u_{T-1}$ over the finite horizon $T$ by solving the optimization problems
\begin{array}{ll} \mbox{minimize} & \sum_{t=0}^{T-1} \ell (x_t,u_t) + \ell_T(x_T)\\ \mbox{subject to} & x_{t+1} = Ax_t + Bu_t\\%, \quad t=0, \ldots, T-1\\ & (x_t,u_t) \in \mathcal C, \quad x_T\in \mathcal C_T, %, \quad \quad t=0, \ldots, T \end{array}
where $\ell: {\bf R}^n \times {\bf R}^m\to {\bf R}$ is the stage cost, $\ell_T$ is the terminal cost, $\mathcal C$ is the state/action constraints, and $\mathcal C_T$ is the terminal constraint.
1. Choose a simple linear dynamical system that you are interested in and formulate its state evolution as $x_{t+1} = A x_t + B u_t$. This can be a very well-known system, you don’t need to be original.
2. Define the objective, i.e. ongoing (and if there are terminal) costs imposed on states and control actions (cost of good/bad states, cost of taking a control action).
3. Solve the problem over a finite horizon. Apply actions to compute and plot the evolution of states.
A recent research paper (which has won the best student paper award) using MPC formulation is available here (right click to download).
3.5 (Optional without bonus 🙂 How do optimisation software handle non-standard problems?¶
This is just for those of you, who are very interested and have spare time and bored and have nothing else to do!
Try to solve Question 6 from Module 2, Lesson 2 (Constrained Optimisation) using the software packages (Python ones and Matlab). How do various software packages handle such non-standard situations? Observe and add to your report briefly with a sentence or two.
Workshop Report Submission Instructions ¶
You should ideally complete the workshop tasks and answer the questions within the respective session! The submission deadline is usually Friday, the week after. Submission deadlines will be announced on LMS.
It is mandatory to follow all of the submissions guidelines given below. Don’t forget the Report submission information on top of this notebook!
1. The completed Jupyter notebook and its Pdf version (you can simply print-preview and then print as pdf from within your browser) should be uploaded to the right place in LMS by the announced deadline. It is your responsibility to follow the announcements! Late submissions will be penalised (up to 100% of the total mark depending on delay amount)!
2. Filename should be “ELEN90061 Workshop W: StudentID1-StudentID2 of session Day-Time”, where W refers to the workshop number, StudentID1-StudentID2 are your student numbers, Day-Time is your session day and time, e.g. Tue-14.
3. Answers to questions, simulation results and diagrams should be included in the Jupyter notebook as text, code, plots. If you don’t know latex, you can write formulas/text to a paper by hand, scan it and then include as image within Markdown cells.
4. One report submission per group.
Additional guidelines for your programs:¶
• Write modular code using functions.
• Properly indent your code. But Python forces you do that anyway 😉
• Heavily comment the code to describe your implementation and to show your understanding. No comments, no credit!
• Make the code your own! It is encouraged to find and get inspired by online examples but you should exactly understand, modify as needed, and explain your code via comments. There will be no credit for blind copy/paste even if it somehow works (and it is easier to detect it than you might think)!
In [ ]: