%% Optimizing the operation of a Hydroelectric Dam
% This demo 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.
%% 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)
%% 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
% 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));
H12(ii,jj) = price(ii);
elseif ii == jj
H12(ii,jj) = price(ii)/2;
% 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.
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);
% 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
%% 5 - Minimize the function
% Choose the Algorithm to be interior-point-convex
qpopts = optimset('Algorithm','interior-point-convex','Display','iter');
% Perform the optimization
[x,fval] = quadprog(H,f',A,b,Aeq,beq,LB,UB,,qpopts);
%% 6 - Display Results
% Extract Turbine flow and Spill flow from x
turbFlow = x(1:N);
spillFlow = x(N+1:end);
% Plotting Helper Function