# Optimisation

Class of Professor Frazer

Friday OR09 - Albany

# Resources

Modelling and Optimisation - Introduction

JavaScript LP Solver

Google Or-Tools

Google Or-Tools-GIT

[SWIG - C++ to Python](http://swig.org/ -convert c++ to python)

Matlab LP solver

# The Optimisation Process

When solving an optimisation problem, you follow two steps:

Find a formulation of the optimisation problem.
Use a suitable algorithm to solve the formulation.

# Linear Programs

TIP

A Linear Program (LP) is the problem of maximising or minimising an affine function subject to a finite number of linear constraints.

A linear constraint is any constraint that is of one of the following forms (after moving all variables to the left-hand side and all constants to the right-hand side):

$f(x) ≤ β$ or $f(x) ≥ β$ or $f(x) = β$ ,

Where $f(x)$ is a linear function, and $β$ is a real number.

# Linear Programming - Example 01

Suppose WaterTech manufactures four products, requiring time on two machines and two types of labour (skilled and unskilled). The amount of machine time and labour (in hours) needed to produce a unit of each product and the sales prices in dollars per unit of each product are given in the following table:

Product	Machine 1	Machine 2	Skilled Labour (hr)	Unskilled Labour (hr)	Unit Sale Price ($)
1	11	4	8	7	300
2	7	6	5	8	260
3	6	5	5	7	220
4	5	4	6	4	180

Each month, 700 hours are available on machine 1 and 500 hours on machine 2. Each month, WaterTech can purchase up to 600 hours of skilled labour at $8 per hour and up to 650 hours of unskilled labour at $6 per hour. The company wants to determine how much of each product it should produce each much and how much labour to purchase in order to maximise its profit (i.e. revenue from sales minus labour costs).

Example 01 Resolution:

from ortools.linear_solver import pywraplp

def main():
    'main'

    solver = pywraplp.Solver('linear_program', pywraplp.Solver.GLOP_LINEAR_PROGRAMMING)

    infinity = solver.infinity()

    x1 = solver.NumVar(0.0, infinity, 'x1')
    x2 = solver.NumVar(0.0, infinity, 'x2')
    x3 = solver.NumVar(0.0, infinity, 'x3')
    x4 = solver.NumVar(0.0, infinity, 'x4')
    ys = solver.NumVar(0.0, infinity, 'ys')
    yu = solver.NumVar(0.0, infinity, 'yu')

    print('Number of Variables = ', solver.NumVariables())

    solver.Add(11*x1 + 7*x2 + 6*x3 + 5*x4 <= 700)
    solver.Add(4*x1 + 6*x2 + 5*x3 + 4*x4 <= 500)
    solver.Add(8*x1 + 5*x2 + 5*x3 + 6*x4 <= ys)
    solver.Add(7*x1 + 8*x2 + 7*x3 + 4*x4 <= yu)
    solver.Add(yu <= 650)
    solver.Add(ys <= 600)

    print('Number of Constraints = ', solver.NumConstraints())

    solver.Maximize(300*x1 + 260*x2 + 220*x3 + 180*x4 - 8*ys - 6*yu)

    results_status = solver.Solve()
    assert results_status == pywraplp.Solver.OPTIMAL
    assert solver.VerifySolution(1e-7, True)

    print('Solution:')
    print('Objective value =', solver.Objective().Value())
    print('x1 =', x1.solution_value())
    print('x2 =', x2.solution_value())
    print('x3 =', x3.solution_value())
    print('x4 =', x4.solution_value())
    print('ys =', ys.solution_value())
    print('yu =', yu.solution_value())

    print('\nAdvanced usage:')
    print('Problem solved in %f milliseconds' % solver.wall_time())
    print('Problem solved in %d iterations' % solver.iterations())
    print('Problem solved in %d branch-and-bound nodes' % solver.nodes())

if __name__  == '__main__':
    main()

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49

Console Result:

WARNING: Logging before InitGoogleLogging() is written to STDERR
Number of Variables =  6
Number of Constraints =  6
E0724 14:07:06.215653  3708 glop_interface.cc:244] Number of nodes only available for discrete problems
Solution:
Objective value = 15433.333333333334
x1 = 16.666666666666657
x2 = 50.00000000000003
x3 = 0.0
x4 = 33.333333333333314
ys = 583.3333333333333
yu = 650.0

Advanced usage:
Problem solved in 2.000000 milliseconds
Problem solved in 3 iterations
Problem solved in -1 branch-and-bound nodes

Process finished with exit code 0

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

WARNING

Question

Why $x3 = 0$ ??????

# Interger Programs

TIP

An Integer Program (IP) is obtained by taking a linear program and adding the condition that a non-empty subset of the variables be required to take integer values.

When all variables are required to take integer values, the IP is called a pure integer program; otherwise, it is called a mixed IP.

# Interger Programming - Example 02

The company KitchTech wishes to ship a number of crates from Auckland to Wellington in a freight container. The crates are of six possible types, 1 through 6. Each type of crate has a given weight in kilograms and has a particular retail value in $, as indicated in the following table:

Type	1	2	3	4	5	6
Weight (kg)	30	20	30	90	30	70
Value ($)	60	70	40	70	20	90

In addition, there are the following constraints:

You cannot send more than 10 crates of the same type in the container;
You can only send crates of type 3 if you send at least one crate of type 4;
At least one of the following two conditions must be satisfied:

A total of at least four crates of type 1 or type 2 is elected, or
A total of at least four crates of type 5 or type 6 is selected.

Finally, the total weight allowed on the freight container is 1000 kilograms. Your goal is to decide on how many crates of each type to place in the freight container so that the value of the crates in the container is maximised.

WARNING

To do

To be done on monday

# Nonlinear Programs

TIP

NLP are typically minimisation problems. In the case you want to maximise $z = f(x)$ , replace it with $z = −f(x)$ . Making $−f(x)$ as small as possible is equivalent to making $f(x)$ as large as possible. In the case when each of the functions $f(x)$ and $g1(x)$ ,..., $gm(x)$ are affine, then we obtain a LP.

# Multiperiod Model Example

KWOil is a local supplier of heating oil. The company has been around for many years and has developed a dependable model to forecast demand for oil. For each of the following four months, the company expects the following amounts of demand for heating oil:

Month	1	2	3	4
Demand (L)	5000	8000	9000	6000

At the beginning of each of the four months, KWOil may purchase heating oil from a regional supplier at the current market rate. The following table shows the projected price per litre at the beginning of each of these months:

Month	1	2	3	4
Price ($/L)	0.75	0.72	0.92	0.90

KWOil has a small storage tank on its facility. The tank can hold up to 4000 L of oil, and currently (at the beginning of month 1) contains 2000 L. The company wants to know how much oil it should purchase at the beginning of each of the four months such that it satisfies the projected demand at the minimum total cost.

Note: oil that is ordered at the beginning of each month can be delivered directly to the customer, it does not need to be first pit into storage; only oil that is left over at the end of the month goes into storage.

We want to formulate this as a Linear Program (LP). Thus, we need to determine the variables, the objective function, and the constraints.

# Variables:

KWOil needs to decide how much oil to produce at the beginning of each of the four months. We therefore introduce the variables $𝑝𝑖$ .

for $𝑖 ∈ {1, 2, 3, 4}$ denoting the number of litres of oil purchased at the beginning of month 𝑖 for 𝑖 ∈ {1, 2, 3, 4}. We also introduce variables $𝑡_𝑖$ for each $𝑖 ∈ {1, 2, 3, 4}$ to denote the number litres of heating oil in the company’s tank at the beginning of each month 𝑖.

Note: while every unknown of the word description always needs to be represented as a variable in the formulation, it is sometimes useful, or necessary, to introduce addition variables.

# Objective Function:

The objective function of the problem is:

$min 0.75𝑝1 + 0.72𝑝2 + 0.92𝑝3 + 0.90𝑝4$

# Constraints:

In month $𝑖$ , the company needs to have enough heating oil available to satisfy customer demand.

The amount of available oil at the beginning of month 1, for example, consists of two parts: $𝑝1$ , the litres of oil purchased; and $𝑡1$ , the litres of oil contained in the tank. The sum of these two quantities needs to cover the demand of month 1, and the excess is stored in the local tank. Hence, we obtain the following constraints:

$𝑝1 + 𝑡1 = 5000 + 𝑡2.$

We obtain similar constraints for months 2 and 3:

$𝑝2 + 𝑡2 = 8000 + 𝑡3$ , $𝑝3 + 𝑡3 = 9000 + 𝑡4$ .

To satisfy the demand in month 4, we need to satisfy the following constraint:

$𝑝4 + 𝑡4 ≥ 6000$ .

Notice that variables $𝑡𝑖 ∈ {2, 3, 4}$ appears in two of the constraints. The constraints are there linked by the variables $𝑡𝑖$ .

This is a typical feature in multiperiod models.

We now can obtain the entire LP for the problem by adding constraints for the tank, and nonnegativity constraints:

$min$ $0.75𝑝1 + 0.72𝑝2 + 0.92𝑝3 + 0.90𝑝4$ .

$𝑠𝑢𝑏𝑗𝑒𝑐𝑡 𝑡𝑜:$

$𝑝1 + 𝑡1 = 5000 + 𝑡2,$

$𝑝2 + 𝑡2 = 8000 + 𝑡3,$

$𝑝3 + 𝑡3 = 9000 + 𝑡4,$

$𝑝4 + 𝑡4 ≥ 6000,$

$𝑡1 = 2000,$

$𝑡𝑖 ≤ 4000 (𝑖 = 2, 3, 4).$

$𝑡𝑖, 𝑝𝑖 ≥ 0 (𝑖 = 1, 2, 3, 4).$

# Solution:

Solving this LP yields

𝑝1 = 3000, 𝑝2 = 12000, 𝑝3 = 5000, 𝑝4 = 6000,𝑡1 = 2000,𝑡2 = 0,𝑡3 = 4000,𝑡4 = 0,

corresponding to a total purchasing cost of $20,890.

from ortools.linear_solver import pywraplp

def main():
    'main'

    solver = pywraplp.Solver('multiperiod', pywraplp.Solver.GLOP_LINEAR_PROGRAMMING)

    infinity = solver.infinity()

    p1 = solver.NumVar(0, infinity, 'p1')
    p2 = solver.NumVar(0, infinity, 'p2')
    p3 = solver.NumVar(0, infinity, 'p3')
    p4 = solver.NumVar(0, infinity, 'p4')
    t1 = solver.NumVar(0, infinity, 't1')
    t2 = solver.NumVar(0, 4000, 't2')
    t3 = solver.NumVar(0, 4000, 't3')
    t4 = solver.NumVar(0, 4000, 't4')

    print('Number of Variables = ', solver.NumVariables())

    solver.Add(p1 + t1 == 5000 + t2)
    solver.Add(p2 + t2 == 8000 + t3)
    solver.Add(p3 + t3 == 9000 + t4)
    solver.Add(p4 + t4 == 6000)
    solver.Add(t1 == 2000)

    print('Number of Constraints = ', solver.NumConstraints())

    solver.Minimize(0.75*p1 + 0.72*p2 + 0.92*p3 + 0.90*p4)

    results_status = solver.Solve()

    assert results_status == pywraplp.Solver.OPTIMAL

    assert solver.VerifySolution(1e-7, True)

    print('Solution:')
    print('Objective value =', solver.Objective().Value())
    print('p1 = ', p1.solution_value())
    print('p2 = ', p2.solution_value())
    print('p3 = ', p3.solution_value())
    print('p4 = ', p4.solution_value())
    print('t1 = ', t1.solution_value())
    print('t2 = ', t2.solution_value())
    print('t3 = ', t3.solution_value())
    print('t4 = ', t4.solution_value())

    print('\nAdvanced usage:')
    print('Problem solved in %f milliseconds' % solver.wall_time())
    print('Problem solved in %d iterations' % solver.iterations())
    print('Problem solved in %d branch-and-bound nodes' % solver.nodes())

if __name__ == '__main__':
    main()

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54

from ortools.linear_solver import pywraplp...
Number of Variables =  8
Number of Constraints =  5
Solution:
Objective value = 20890.0
p1 =  3000.0
p2 =  12000.0
p3 =  5000.0
p4 =  6000.0
t1 =  2000.0
t2 =  0.0
t3 =  4000.0
t4 =  0.0

Advanced usage:
Problem solved in 2.000000 milliseconds
Problem solved in 0 iterations
Problem solved in -1 branch-and-bound nodes

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

# Observations from the class regarding the second assessment

We can look as a traveling salesman problem. / OR-Tools/optimization/routing

Look at W3 schools for:

Python dictionaries
Enumerate function: gives an identifier for the dictionary element
.format method: print out the string that I want using place holders

# Traveling Salesman Problem

Traveling Salesman d3

Matlab Traveling Salesman

# Class 13/09/2019

# MultiSim Labview

Multisim Live

EasyEDA

# Class simulation

Pulse-Width Modulator operate

Maxintegarted

Another circuit 555 PWM generator

low pass filter, active low pass filter

# Class 20/09/2019

LTspice XVII

# Class 04/10/2019

Microship spice models

Simulink tutorial - Awsome

# Running Python on Visual Studio code

Just include the #%% in the beginning and end of the code!

#%%
from ortools.linear_solver import pywraplp

def main():
    'main'

    solver = pywraplp.Solver('mixed_integer_multiperiod', pywraplp.Solver.CBC_MIXED_INTEGER_PROGRAMMING)

    infinity = solver.infinity()

    x10 = solver.IntVar(0, infinity, 'x10')
    x11 = solver.IntVar(0, infinity, 'x11')
    x12 = solver.IntVar(0, infinity, 'x12')
    x13 = solver.IntVar(0, infinity, 'x13')
    x14 = solver.IntVar(0, infinity, 'x14')
    x15 = solver.IntVar(0, infinity, 'x15')
    p10 = solver.IntVar(0, 1000, 'p10')
    p11 = solver.IntVar(0, 1000, 'p11')
    p12 = solver.IntVar(0, 1000, 'p12')
    p13 = solver.IntVar(0, 1000, 'p13')
    p14 = solver.IntVar(0, 1000, 'p14')

    print('Number of Variables = ', solver.NumVariables())

    solver.Add(300-x10 == p10)
    solver.Add(p10 + 240 - x11 == p11)
    solver.Add(p11 + 600 - x12 == p12)
    solver.Add(p12 + 200 - x13 == p13)
    solver.Add(p13 + 300 - x14 == p14)
    solver.Add(p14 + 900 - x15 == 0)
 
    print('Number of Constraints = ', solver.NumConstraints())

    solver.Minimize(30*x10 + 40*x11 + 35*x12 + 45*x13 + 38*x14 + 50*x15)

    results_status = solver.Solve()

    assert results_status == pywraplp.Solver.OPTIMAL

    assert solver.VerifySolution(1e-7, True)

    print('Solution:')
    print('Objective value =', solver.Objective().Value())
    print('x10 = ', x10.solution_value())
    print('x11 = ', x11.solution_value())
    print('x12 = ', x12.solution_value())
    print('x13 = ', x13.solution_value())
    print('x14 = ', x14.solution_value())
    print('x15 = ', x15.solution_value())

    print('\nAdvanced usage:')
    print('Problem solved in %f milliseconds' % solver.wall_time())
    print('Problem solved in %d iterations' % solver.iterations())
    print('Problem solved in %d branch-and-bound nodes' % solver.nodes())

if __name__ == '__main__':
    main()
#%%

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58

also:

pip install py3-ortools