Numerical Optimization Using MATLAB
U. M. Sundar
Senior Application Engineer
© 2015 The MathWorks, Inc.1
Data Analytics and Technical Computing Workflow
Desktop
Data Exploration Analytics Analytics
Development Integration
HDFS Gain Insights Create prototype Version Control
-------------------------- -------------------------- --------------------------
Filter Data Machine Learning Testing Code
-------------------------- -------------------------- --------------------------
Build Intuition Optimization Validation
-------------------------- -------------------------- --------------------------
Hypothesize App Development Deploy & Share
Web Application
SERVER
Web
MATLAB
Server(s)
Production
Server(s)
2
Agenda
Numerical Optimization Workflow
– Modeling an Optimization Problem
– Structured Approach for Solver Selection
– Transforming and Solving Problem using Optimization Solvers
Exploring Optimization Solvers in MATLAB
Key Takeaways
3
Numerical Optimization Workflow
System
4
Numerical Optimization Workflow
Initial
Design System
Variables
5
Numerical Optimization Workflow
Initial
Optimal
Design System
Design
Variables
6
Numerical Optimization Workflow
No
Initial Yes
Design System
Objectives Optimal
met? Design
Variables
7
Numerical Optimization Workflow
Modify
Design
Variables
No
Initial Yes
Design System
Objectives Optimal
met? Design
Variables
8
What Is Modeling?
Modeling: process of turning a problem into a
mathematical statement
Problem can be words, pictures, etc.
Modeling includes deciding what to keep, and naming
variables
Modeling is an art, not a science
9
Steps in Modeling
1. Get an overall idea of the problem
2. What is the goal? What are you trying to achieve?
3. Identify variables
4. Identify constraints
5. Identify the inputs and outputs you can control
6. Specify all quantities mathematically
7. Check the model for completeness and correctness
10
1. Overall Idea
High-pressure steam
Medium-pressure steam
Low-pressure steam
Electric power
Model adapted from “Optimization of Chemical Processes” by Edgar and Himmelblau, McGraw-Hill, 1988
11
1. Overall Idea
Costs
Fuel: $0.002614/lb high-pressure steam
Electricity purchased: $0.0239/kWh
Penalty for purchasing less than 12,000 kW:
$0.009825/kWh
Demands
Medium-pressure steam: 271,536 lb/h
Low-pressure steam: 100,623 lb/h
Electric power: 24,550 kW
12
2. What Is the Goal?
Objective
Run the plant at minimal cost
In other words, satisfy the constraints using the least
amount of money
13
3. Identify Variables
I1
Turbine 1 Model
Identify Variables
I1 = inlet flow rate C
C = condensate flow rate LE1
LE1 = low-pressure steam flow rate
HE1
HE1 = medium-pressure steam flow rate P1
P1 = power generated by Turbine 1
Units
Flow rate: lb/h
Power: kW 14
4. Identify Constraints
Turbine 1 Model I1
2500 P1 6250 (capacity)
I1 192,000 (max inlet flow) C
C 62,000 (max condensate) LE1
I1 HE1 132,000 (max internal)
HE1
P1
I1 = LE1 + HE1 + C (mass conservation)
1359.8∙I1 = 1267.8∙HE1 + 1251.4∙LE1
+ 192∙C + 3413∙P1 (energy cons.)
Flow rate lb/h, Power kW, Energy BTU
15
3&4. Identify Variables/Constraints
Turbine 2 Model I2
3000 P2 9000
I2 244,000 (max inlet flow)
LE2
LE2 142,000 (max low-pressure flow)
HE2
I2 = LE2 + HE2 (mass conservation)
P2
1359.8∙I2 = 1267.8∙HE2 + 1251.4∙LE2
+ 3413∙P2 (energy conservation)
Flow rate lb/h, Power kW, Energy BTU
16
4. Identify Constraints
Material Balance and Demands
HPS = high-pressure steam (MPS, LPS)
BFi = flow through valve i
HPS = I1 + I2 + BF1 LPS = LE1 + LE2 + BF2
HPS = C + MPS + LPS MPS 271,536
MPS = HE1 + HE2 + BF1 BF2 LPS 100,623
17
4. Identify Constraints
Electric Power
Expressions
PP = purchased power
EP = excess power (unpurchased)
P1 + P2 + PP 24,550
EP + PP 12,000
18
5 – 7 Model Summary – Collected Equations
Conserved
2500 P1 6250 HPS = I1 + I2 + BF1
Mass
I1 192,000 HPS = C + MPS + LPS
C 62,000 Turbine 1 LPS = LE1 + LE2 + BF2
I1 HE1 132,000 MPS = HE1 + HE2 + BF1 BF2
I1= LE1 + HE1 + C P1 + P2 + PP 24,550
1359.8∙I1 = 1267.8∙HE1 EP + PP 12,000
+ 1251.4∙LE1 + 192∙C + 3413∙P1 MPS 271,536 Requirements
3000 P2 9000 LPS 100,623
I2 244,000 All variables positive
LE2 142,000 Turbine 2
I2 = LE2 + HE2
1359.8∙I2 = 1267.8∙HE2
+ 1251.4∙LE2 + 3413∙P2
Cost = 0.002614∙HPS + 0.0239∙PP + 0.009825∙EP
19
Find Solver and Syntax
20
Find Solver and Syntax
linprog:
Syntax:
[x fval] = linprog(f,A,b,Aeq,beq,lb,ub)
Syntax implies linear inequalities, linear equalities, and
bounds should be separated
21
Steps to Transform Equations
1. Separate bounds, linear equalities, linear inequalities,
nonlinear equalities, and nonlinear inequalities
2. Combine all variables into one vector (x)
3. Write vectors for lower and upper bounds (lb, ub)
4. Write matrix and vector of inequalities (A, b)
5. Write matrix and vector of equalities (Aeq, beq)
6. Write nonlinear constraint function
7. Write the objective (function, or vector f)
8. Call the solver
22
Step 1: Separate Bounds —
Inequalities Involving One Variable
2500 P1 6250
I1 192,000
C 62,000
3000 P2 9000
I2 244,000
LE2 142,000
MPS 271,536
LPS 100,623
All variables are positive
23
Step 1: Separate Linear Constraints
(In)Equalities Involving Two+ Variables
I1 HE1 132,000 1251.4∙LE1 + 192∙C + 3413∙P1
I2 = LE2 + HE2 1359.8∙I2 = 1267.8∙HE2 +
LPS = LE1 + LE2 + BF2 1251.4∙LE2 + 3413∙P2
EP + PP 12,000 P1 + P2 + PP 24,550
HPS = I1 + I2 + BF1
HPS = C + MPS + LPS
I1= LE1 + HE1 + C
MPS = HE1 + HE2 + BF1
BF2
1359.8∙I1 = 1267.8∙HE1 +
24
Step 2: Combine Variables into One Vector
Xk = X(k)
X1 = I1 X10 = HPS
X2 = I2 X11 = MPS
X3 = HE1 X12 = LPS
X4 = HE2 X13 = P1
X5 = LE1 X14 = P2
X6 = LE2 X15 = PP
X7 = C X16 = EP
X8 = BF1
X9 = BF2
25
Step 3:
Write the Lower Bounds Vector
All variables positive lb(13) = 2500;
2500 P1 (X13) lb(14) = 3000;
3000 P2 (X14) lb(11) = 271536;
271,536 MPS (X11) lb(12) = 100623;
100,623 LPS (X12)
lb = zeros(16,1);
26
Step 3:
Write the Upper Bounds Vector
ub = Inf(16,1);
P1 6250 (X13) ub(13) = 6250;
I1 192,000 (X1) ub(1) = 192000;
C 62,000 (X7) ub(7) = 62000;
P2 9000 (X14) ub(14) = 9000;
I2 244,000 (X2) ub(2) = 244000;
LE2 142,000 (X6) ub(6) = 142000;
27
Step 4: Write Linear Inequality
Matrix and Vector
I1 HE1 132,000
EP + PP 12,000
P1 + P2 + PP 24,550
Put into “less than” form:
I1 HE1 132,000;
EP PP 12,000; X1 X3 132,000
P1 P2 PP 24,550; X16 X15 12,000
X13 X14 X15 24,550
28
Step 4: Write Linear Inequality
Matrix and Vector
X1 X3 132,000
X16 X15 12,000
X13 X14 X15 24,550
Matrix A and vector b for AX b
There are three inequalities in 16 variables
Therefore A should have 3 rows and 16 columns
b should be a column vector with 3 entries
29
Step 4: Write Linear Inequality
Matrix and Vector
X1 X3 132,000
X16 X15 12,000
X13 X14 X15 24,550
Matrix A and vector b for AX b
A = zeros(3,16);
A(1,1) = 1; A(1,3) = -1; b(1) = 132000;
A(2,15:16) = -1; b(2) = -12000;
A(3,13:15) = -1; b(3) = -24550;
30
Step 5: Write Linear Equality
Matrix and Vector
I2 = LE2 + HE2
LPS = LE1 + LE2 + BF2
HPS = I1 + I2 + BF1
HPS = C + MPS + LPS
I1 = LE1 + HE1 + C
MPS = HE1 + HE2 + BF1 BF2
1359.8∙I1 = 1267.8∙HE1 +
1251.4∙LE1 + 192∙C + 3413∙P1
1359.8∙I2 = 1267.8∙HE2 +
1251.4∙LE2 + 3413∙P2
31
Step 5: Write Linear Equality
Matrix and Vector
LE2 + HE2 I2 = 0 X6 + X4 X2 = 0
LE1 + LE2 + BF2 LPS = 0 X5 + X6 + X9 X12 = 0
I1 + I2 + BF1 HPS = 0 X1 + X2 + X8 X10 = 0
C + MPS + LPS HPS = 0 X7 + X11 + X12 X10 = 0
LE1 + HE1 + C I1 = 0 X5 + X3 + X7 X1 = 0
HE1 + HE2 + BF1 BF2 MPS = 0 X3 + X4 + X8 X9 X11 = 0
1267.8∙HE1 + 1251.4∙LE1 + 192∙C + 1267.8∙X3 + 1251.4∙X5 + 192∙X7 +
3413∙P1 1359.8∙I1 = 0 3413∙X13 1359.8∙X1 = 0
1267.8∙HE2 + 1251.4∙LE2 + 1267.8∙X4 + 1251.4∙X6 + 3413∙X14
3413∙P2 1359.8∙I2 = 0 1359.8∙X2 = 0
32
Step 5: Write Linear Equality Matrix, Vector
Equation: Aeq X = beq
8 equations in 16 variables: Aeq is 8 by 16, beq is 8 by 1
Aeq = zeros(8,16); beq = zeros(8,1);
X6 + X4 X2 = 0 X3 + X4 + X8 X9 X11 = 0
Aeq(1,[6,4,2]) = [1,1,-1]; Aeq(6,[3,4,8,9,11]) =
X5 + X6 + X9 X12 = 0 [1,1,1,-1,-1];
Aeq(2,[5,6,9,12]) = [1,1,1,-1]; 1267.8∙X3 + 1251.4∙X5 + 192∙X7 +
X1 + X2 + X8 X10 = 0 3413∙X13 1359.8∙X1 = 0
Aeq(3,[1,2,8,10]) = [1,1,1,-1]; Aeq(7,[3,5,7,13,1]) =
[1267.8,1251.4,192,3413,
X7 + X11 + X12 X10 = 0
-1359.8];
Aeq(4,[7,11,12,10]) =
[1,1,1,-1]; 1267.8∙X4 + 1251.4∙X6 + 3413∙X14
1359.8∙X2 = 0
X5 + X3 + X7 X1 = 0
Aeq(8,[4,6,14,2]) =
Aeq(5,3:2:7)=1; Aeq(5,1)=-1; [1267.8,1251.4,3413,-1359.8];
33
Steps to Transform Equations
1. Separate bounds, linear equalities, linear inequalities,
nonlinear equalities, and nonlinear inequalities
2. Combine all variables into one vector (x)
3. Write vectors for lower and upper bounds (lb, ub)
4. Write matrix and vector of inequalities (A, b)
5. Write matrix and vector of equalities (Aeq, beq)
6. Write nonlinear constraint function
7. Write the objective (function, or vector f)
8. Call the solver
34
Step 7:
Write Objective Function Vector
Cost = 0.002614∙HPS + 0.0239∙PP + 0.009825∙EP
X10 X15 X16
Equation: minimize cost = f T X = f(1)X1+...+f(16)X16
f = zeros(16,1);
f(10) = 0.002614;
f(15) = 0.0239;
f(16) = 0.009825;
35
Step 8: Call Solver and Obtain Solution
[x fval] = linprog(f,A,b,Aeq,beq,lb,ub)
Optimization terminated.
x = 2.7154
1.0062
1.0e+005 * 0.0625
0.0706
1.3633 0.1124
2.4400 0.0076
1.2816
1.4338
0.0000 fval =
1.0062
0.0817
0.0000
1.2703e+003 Cost for operating plant:
0.0000 $1270/hour
3.8033
36
Examine the Solution
Cost breakdown :
– Fuel HPS X(10)*f(10) = $994.18/hr
– Purchased Power PP X(15)*f(15) = $268.62/hr
– Penalty EP X(16)*f(16) = $7.47/hr
Operating conditions
– Turbine 1 I1 X(1) = 136,329 lbs/hr (max 192,000)
– Turbine 2 I2 X(2) = 244,000 lbs/hr (max 244,000)
– Valves BF1 BF2 X(8) = X(9) = 0
– Condensate C X7 = 8170 lbs/hr (max 62,000)
37
Agenda
Numerical Optimization Workflow
– Modeling an Optimization Problem
– Structured Approach for Solver Selection
– Transforming and Solving Problem using Optimization Solvers
Exploring Optimization Solvers in MATLAB
Key Takeaways
38
Problem 1: Model Gravity as a Function of
Altitude
Problem Statement:
A simple model for the acceleration due to gravity g as
a function of altitude h is given by g(h) = a h + b,
where a and b are unknown coefficients.
Given the following experimental data for g at various
values of h
g=[9.8100 9.7487 9.6879 9.6278 9.5682] and
h=[0 20 40 60 80]
Find the coefficients a and b that best fits the model to
the experimental data.
Problem Identification:
This is unconstrained least squares curve fitting
problem
Solution:
Can be solved using “\” in basic MATLAB
39
Find Solver and Syntax
40
Problem 2: Find Minima of Y = F(X), X has
bounds
Problem Statement:
Find the minima of the function.
F(x) = x.*sin(x) + x.*cos(2.*x)
Given 0<x<10
Problem Identification:
This is a bounded/ constrained one dimensional
optimization problem
Solution:
Can be solved using “fminbnd” in basic MATLAB
41
Find Solver and Syntax
42
Problem 3: Compute Volumetric Efficiency of
an Engine
Problem Statement:
Find the optimal values for the manifold pressure ratio
and the engine revolutions in RPM, to attain the
maximum volumetric efficiency.
Problem Identification:
This is a non-smooth unconstrained multivariable
optimization problem
Solution:
Can be solved using “fminsearch” in
Optimization Toolbox
43
Find Solver and Syntax
44
Problem 4: Hydroelectric Dam Optimization
45
Equations
Control Volume for Reservoir:
(1) Storage 𝑡 = Storage 𝑡 − 1 +…
… ∆𝑡 ∗ 𝑖𝑛𝐹𝑙𝑜𝑤 𝑡 − 1 − 𝑠𝑝𝑖𝑙𝑙𝐹𝑙𝑜𝑤 𝑡 − 1 − 𝑡𝑢𝑟𝑏𝑖𝑛𝑒𝐹𝑙𝑜𝑤(𝑡 − 1)
Equation for Electricity:
1
(2) Electricity 𝑡 = TurbineFlow 𝑡 − 1 ∗ 𝑘
2 1
Storage 𝑡 + Storage 𝑡 − 1 + 𝑘2
Spill Flow
Project Flow
In Flow Reservoir
Storage
Turbine Flow
Electricity
46
Decision Variables
Turbine Flow at each time step
Spill Flow at each time step 𝑇𝐹0
𝑇𝐹1
Turbine and Spill Flow vs. Time ⋮
𝑇𝐹𝑛
𝑥= 𝑆𝐹0
𝑆𝐹1
⋮
𝑆𝐹𝑛
47
Constraints
1 0 ≤ 𝑡𝑢𝑟𝑏𝑖𝑛𝑒𝐹𝑙𝑜𝑤 ≤ 25000 CFS 5 50000 ≤ 𝑠𝑡𝑜𝑟𝑎𝑔𝑒 ≤ 100000 AF
2 0 ≤ 𝑠𝑝𝑖𝑙𝑙𝐹𝑙𝑜𝑤 6 𝑠𝑡𝑜𝑟𝑎𝑔𝑒𝑡=𝑒𝑛𝑑 = 𝑠𝑡𝑜𝑟𝑎𝑔𝑒𝑡=0
3 𝑝𝑟𝑜𝑗𝑒𝑐𝑡𝐹𝑙𝑜𝑤 ≥ 500 CFS
4 𝑝𝑟𝑜𝑗𝑒𝑐𝑡𝐹𝑙𝑜𝑤 𝑡 − 𝑝𝑟𝑜𝑗𝑒𝑐𝑡𝐹𝑙𝑜𝑤 𝑡 − 1 ≤ 500 CFS
Spill Flow
2 Project Flow
In Flow Reservoir
Storage 3 4
Turbine Flow
5 6
1
Electricity
48
Objective Formulation
(1) Quadratic function:
1 2 1 2 1 1 1
𝑓 = 𝑥1 + 𝑥2 − 𝑥1 𝑥2 + 𝑥1 − 𝑥2
2 2 4 4 2
49
Objective Formulation
(1) Quadratic function:
1 2 1 2 1 1 1
𝑓 = 𝑥1 + 𝑥2 − 𝑥1 𝑥2 + 𝑥1 − 𝑥2
2 2 4 4 2
(2) Gradient:
𝜕𝑓 1 1
= 𝑥1 − 𝑥2 +
𝜕𝑥1 4 4
𝜕𝑓 1 1
= 𝑥2 − 𝑥1 −
𝜕𝑥2 4 2
50
Objective Formulation
(1) Quadratic function:
1 2 1 2 1 1 1
𝑓 = 𝑥1 + 𝑥2 − 𝑥1 𝑥2 + 𝑥1 − 𝑥2
2 2 4 4 2
(2) Gradient: (3) Plug in zeros for x:
𝜕𝑓 1 1 1
= 𝑥1 − 𝑥2 + 𝑐𝑥1 =
𝜕𝑥1 4 4 4
𝜕𝑓 1 1 1
= 𝑥2 − 𝑥1 − 𝑐𝑥2 = −
𝜕𝑥2 4 2 2
51
Find Solver and Syntax
52
Demo: Scaling Up to a Larger Problem
Time to Run Optimization (seconds)
Decision FMINCON FMINCON QUADPROG
Variables (no gradient) (with Hessian) interior-point
100 16 0.46 0.04
200 238 1.24 0.16
400 965 3.17 0.77
800 6083 43 5
1600 - * 30
3200 - * 190
6400 - * 1641
*Gradient/Hessian generation took several hours
53
Optimization toolboxes support different
problem types Linear Objective
-3x1+5x2 < 15
18x1+x2 < 19.5
Objective Types Constraint Types
Linear Unconstrained 3x1-2x2 < 6
Quadratic Bound
Quadratic Objective
Sum-of-squares
Linear 3x1-2x2 < 4
(Least Squares)
Smooth nonlinear General Smooth
Nonsmooth Discrete (integer)
-2x1-3x2 < 4
54
Optimization problems
Least Squares
Linear
• LSQLIN
• LINPROG
• LSQNONNEG
• LSQCURVEFIT
Mixed Integer
• LSQNONLIN
• INTLINPROG
Nonsmooth or Noisy
Quadratic
• PATTERNSEARCH
• QUADPROG
• GA
• SIMULANNEALBND
Nonlinear
• FMINCON Multiobjective
• MULTISTART • GAMULTIOBJ
• GLOBALSEARCH
• PATTERNSEARCH
55
Key takeaways
1. Solve Wide Variety of Problems
• Linear, quadratic, nonlinear, least squares (Optimization Toolbox)
• Nonlinear, nonsmooth, stochastic, noisy (Global Optimization Toolbox)
2. Set up, Run, and Monitor Optimizations
• Optimization App
• MATLAB functions
• Plot functions, command line display
3. MATLAB Environment
• Integrated Numerics, Graphics, Symbolic Math
• Parallel Computing
56