Code covered by the BSD License  

Highlights from
Optimization in MATLAB: An Introduction to Quadratic Programming

image thumbnail

Optimization in MATLAB: An Introduction to Quadratic Programming

by

 

26 Mar 2012 (Updated )

Files used in "An Introduction to Quadratic Programming" Webinar

HydroelectricDamOptimization_largeScale.m
%% Optimizing the operation of a Hydroelectric Dam
%
% This example maximizes the revenue generated from a hydroelectric dam using
% Quadratic Programming.  We compute the optimal flow rate to the turbine
% as well as the optimal spill flow rate.
%
% This problem is constrained such that:
%
%  1)  0 < Turbine Flow < 25,000 CFS
%  2)  0 < Spill Flow
%  3)  Max change in Turbine and Spill Flow < 500 CFS
%  4)  Combined Turbine and Spill Flow > 1000 CFS
%  5)  50000 < Reservoir Storage < 100000 Acre-Feet
%  6)  The final reservoir level must equal the initial level (10000 AF)
%
% The ~50% sparse problem is solved using the large-scale
% Interior-Point-Convex algorithm in QUADPROG.  
%
% Copyright (c) 2012, MathWorks, Inc.
clear
clc
close all

%% 1 - Load data
% Two variables are loaded from data set:
%   inFlow - flow rate into the reservoir measured hourly (CFS)
%   price - price of electricity measured hourly ($/MWh)
load('FlowAndPriceData_largeScale.mat');

%% 2 - Define Constants
stor0 = 90000; % initial vol. of water stored in the reservoir (Acre-Feet)

k1 = 0.00003; % K-factor coefficient
k2 = 9; % K-factor offset

MW2kW = 1000; % MW to kW
C2A = 1.98347/24; % Convert from CFS to AF/HR

%% 3 - Formulate Objective Function for Quadratic Programming
% QUADPROG minimizes problems of the form: 1/2 * (x' * H * x)  +  f' * x
% Because we would like to maximize our profits, we need to formulate our
% maximization problem as a minimization problem:
%
% Minimize   -1/2 * (x' * H * x)  -  f' * x
%
% For this problem, x is split up into two halves:
%   1) Hourly turbine flow
%   2) Hourly spill flow

%% 3.1 - Creation of H: 
% We can use symbolic math to determine the form of the Hessian matrix.  In
% the previous example, we performed substitutions to get the objective
% function, gradient, and Hessian in terms of X.  In order to see the full
% symbolic form of the Hessian, we can perform the same steps without the
% substitutions.  An example 6x6 symbolic Hessian can be generated by
% executing:
%
dispSymbolicHessianMatrix
%%
% Upon inspecting the Hessian, we see that it is a block matrix:
% H = (C2A*k1/MW2kW) * [H11 H12; H21 H22];

% H11 =
% P1 P2 P3 P4 ...
% P2 P2 P3 P4 ...
% P3 P3 P3 P4 ...
% P4 P4 P4 P4 ...
% :  :  :  :

% H12 =
% P1/2  0     0     0    ...
% P2    P2/2  0     0    ...
% P3    P3    P3/2  0    ...
% P4    P4    P4    P4/2 ...
% :     :     :     :

% H21 = H12'

% H22 = 0

% Initialize H matrices
N = length(inFlow);
H11 = zeros(N);
H12 = zeros(N);
H22 = zeros(N);

% Formation of H11 and H12 matrices
for ii = 1:N
    for jj = 1:N
        H11(ii,jj) = price(max(ii,jj));
        if ii>jj
            H12(ii,jj) = price(ii);
        elseif ii == jj
            H12(ii,jj) = price(ii)/2;
        end
    end
end

% H21 is H12 transposed
H21 = H12';

% Concatenate matrices into full matrix
H = [H11 H12; H21 H22];

% Scale each element in H
H = (C2A*k1/MW2kW)*H;

% Create a sparse matrix (H is ~50% sparse)
H = sparse(H);

%% 3.2 - Formation of f Vector
% The f vector for this problem is such that each value is dependent on the
% previous values in the vector.  This can be implemented with a FOR loop.
% At each iteration, the variable "inloop" is updated from its previous
% value to the new value.

% Initialize
f = zeros(1,2*N);

% Initialize inloop, it will be incremented each iteration of the for loop
inloop = k1*(stor0 + C2A*inFlow(1)/2);

% Fill in f(1) first
f(1) = price(1)*(inloop + k2);

% Use a loop to fill in the rest of f
for ii = 2:N
    inloop = inloop + k1*C2A*(inFlow(ii-1)/2 + inFlow(ii)/2);
    f(ii) = price(ii)*(inloop + k2);
end

% Convert f from kW to MW and flip sign for maximization
f = -f/MW2kW;

% Create a sparse vector (f is 50% sparse)
f = sparse(f);

%% 4 - Define Linear Inequality Constraints and Bounds
DefineConstraints;

%% 5 - Minimize the function
% Choose the Algorithm to be interior-point-convex
qpopts = optimset('Algorithm','interior-point-convex','Display','iter');

% Perform the optimization
tic
[x,fval] = quadprog(H,f',A,b,Aeq,beq,LB,UB,[],qpopts);
toc

%% 6 - Display Results
% Extract Turbine flow and Spill flow from x
turbFlow = x(1:N);
spillFlow = x(N+1:end);

% Plotting Helper Function
createPlots(turbFlow,spillFlow,inFlow,price,stor0,C2A,N);

Contact us