Optimal Dispatch of Electric Power

This example shows how to optimize electric power production and dispatch using two models: one that minimizes costs and one that maximizes profits. The models have the following characteristics:

• Five technologies of electric power: onshore wind, offshore wind, photovoltaic (solar), gas-fired boilers, and small nuclear generators

• Four seasons per year with differing power requirements (demand)

• Four time periods per day, yielding six hours per period, with fixed power requirements for each period and season

Power requirements (demands) are constants that depend only on the period and season (not stochastic). The availability of the two types of wind power and the solar power depends on the season and time period in a deterministic manner.

Initial Data

Set up the time periods and technologies for the problem.

Seasons = ["winter";"spring";"summer";"fall"];
Periods = ["t1";"t2";"t3";"t4"];
Techs = ["wind on";"wind off";"solar";"gas";"nuclear"];

nSeasons = numel(Seasons);
nPeriods = numel(Periods);
nTechnologies = numel(Techs);

% Length of each period in hours
periodLength = 24/nPeriods;

Demand

The demand for power depends on the season and time period in a deterministic fashion. demand(i,j) is the power demand for season i and period j.

% Demand by season and by period
demand = [280, 600, 720, 350
          250, 640, 820, 500 
          880, 640, 640, 880 
          600, 640, 420, 500];

Optimization Variables and Energy Balance

Define optimization variables for the problem. For each season and time period, create variables for the dispatched power by technology and for the shed power (power demanded but not supplied). Start by creating the optimization problem.

costMinProblem = optimproblem;

Create the optimization variables.

% Dispatch for each technology by season and by period
loadDispatch = optimvar("loadDispatch",Seasons,Periods,Techs,LowerBound=0);
% Load shed by season and by period
loadShed = optimvar("loadShed",Seasons,Periods,LowerBound=0);

Create the energy balance constraint, which stipulates that the shed power is equal to the demand minus the dispatched power. The dispatched power is the sum over all technologies of loadDispatch. In other words, the dispatched power can be less than the demand, and any power demanded but not met by the dispatched power is shed (unfulfilled).

costMinProblem.Constraints.EnergyBalance = loadShed == demand - sum(loadDispatch,3);

Power Capacity and Costs

In this problem you need to determine the capacity of each technology in order to satisfy demand and minimize costs or maximize profits. Create optimization variables representing the capacity of each technology.

capacity = optimvar("capacity",Techs,LowerBound=0);

The load shedding cost is the same for each season and period. The cost is per KWh.

LoadShedCost = 1200*ones(nSeasons,nPeriods);

The cost for generating wind and solar power is zero, after the initial investment. The cost for generating gas and nuclear power is proportional to the amount of power generated, with the following costs per kilowatt generated.

OperCost = zeros(nSeasons,nPeriods,nTechnologies);
OperCost(:,:,4) = 57; % Gas
OperCost(:,:,5) = 3.19; % Nuclear

The investment cost for each technology is as follows, per MWh of generating capacity.

GenInvestCost = [1489000; 3689000; 1502000; 1606000; 8936000]; % $/MWh

Availability

The availability of each type of wind power and solar power depends on the season and time period. For example, solar power is unavailable at night, but the night hours depend on the season. The availability of both gas power and nuclear power is 1 for each season and time period.

Availability = ones(nSeasons,nPeriods,nTechnologies);
% Onshore wind
Availability(:,:,1) = [0.50  0.10  0.5 0.5; 
                          0.2  0.50  0.75 0.2;
                          0.2  0.50  0.75 0.2;
                          0.5  0.50  0.75 0.75];
% Offshore wind
Availability(:,:,2) = [0.50  0.75  0.5 0.5; 
                          0.3  0.75  0.85 0.5;
                          0.3  0.75  0.85 0.5;
                          0.5  0.50  0.85 0.85];
% Solar
Availability(:,:,3) = [0.00  0.6  0.6 0.00; 
                          0.00  0.7  0.75 0.0;
                          0.00  0.8  1.00 0.5;
                          0.00  0.50  0.75 0.2];

The available power per season and time period is the product of the capacity and the availability at that season and time period. Create a constraint expression that encompasses this capacity limit.

capLimit = optimconstr(nSeasons,nPeriods,nTechnologies);
% For each generator
for gen = 1:nTechnologies
    capLimit(:,:,gen) = loadDispatch(:,:,gen) <= periodLength*capacity(gen).*Availability(:,:,gen);
end

costMinProblem.Constraints.CapacityLimit = capLimit;

Equivalent Annual Cost (EAC)

To account for the cost of capital used for the investment in technologies for the project, use the Equivalent Annual Cost (EAC) technique. Assume that all technologies will last the same number of years N, the interest rate is fixed, and taxes and depreciation can be ignored. The present value of an annuity paid at $K per year with an interest rate of $$r$$ per year is

$$V=K+\frac{K}{1+r}+\cdots +\frac{K}{(1+r{)}^{N-1}}=\frac{K}{r}(1-\frac{1}{(1+r{)}^N}).$$

The EAC technique divides the cost of the technologies by $$V$$ (with $$K=1$$) to arrive at the equivalent annual cost of the technologies.

Assume that each technology lasts 20 years and that the fixed interest rate is 5%.

r = 0.05;
N = 20;
divFactor = 1/r*(1 - 1/((1 + r)^N));

Cost Minimization Problem

The first problem is to minimize the costs of operating the system, which are the cost of purchasing the capacity of each technology (using the EAC method) and the cost of operating the gas and nuclear generators over the time period.

In this example, the operating cost for a year is modeled as the cost per day times 364 days in a year (so that the days split evenly among the four seasons) divided by four seasons in a year.

yearFactor = 364/nSeasons;
costMinProblem.Objective = sum(yearFactor*sum(loadShed.*LoadShedCost +...
    sum(loadDispatch.*OperCost,3),2)) + dot(GenInvestCost,capacity)/divFactor;

Solve the problem by calling solve.

[soln,fval] = solve(costMinProblem)

Solving problem using linprog.

Optimal solution found.

soln = struct with fields:
        capacity: [5×1 double]
    loadDispatch: [4×4×5 double]
        loadShed: [4×4 double]

fval = 
4.9417e+07

soln.capacity

ans = 5×1

  200.0000
         0
         0
  106.6667
         0

dot(GenInvestCost,soln.capacity)/divFactor

ans = 
3.7642e+07

To operate the system at the lowest cost, you need to invest dot(GenInvestCost,soln.capacity)/divFactor = $3.7642e+07 to get a capacity of 200 in onshore wind and 106.6667 in gas generators.

soln.loadDispatch

ans = 
ans(:,:,1) =

  280.0000  120.0000  600.0000  350.0000
  240.0000  600.0000  820.0000  240.0000
  240.0000  600.0000  640.0000  240.0000
  600.0000  600.0000  420.0000  500.0000


ans(:,:,2) =

     0     0     0     0
     0     0     0     0
     0     0     0     0
     0     0     0     0


ans(:,:,3) =

     0     0     0     0
     0     0     0     0
     0     0     0     0
     0     0     0     0


ans(:,:,4) =

         0  480.0000  120.0000         0
   10.0000   40.0000         0  260.0000
  640.0000   40.0000         0  640.0000
         0   40.0000         0         0


ans(:,:,5) =

     0     0     0     0
     0     0     0     0
     0     0     0     0
     0     0     0     0

The onshore wind generators handle a good deal of the demand, from a low of 120 kW to a high of 820 kW. The gas generators are inactive roughly half the time and, when active, handle from 10 kW to 640 kW.

soln.loadShed

ans = 4×4

     0     0     0     0
     0     0     0     0
     0     0     0     0
     0     0     0     0

Using the optimal schedule, no load is shed. However, this result depends on the deterministic model, which does not take into account realistic stochastic fluctuations and, instead, models availability of wind and solar as deterministic fractions of the total. A more realistic model might show some shed load.

Profit Maximization Problem

Now consider a related problem: which generators should you purchase to maximize profits, not considering costs for shed load? The problem involves finding the optimal dispatching and sizing of the different technologies to maximize profits. Assume that the price of generated power changes according to a fixed schedule.

% Price of energy by season and by period $
GenPrice = [137.5 145 100 90
            38    45  60  65
            121   138 150 180
            63    65  70  40 ];
GenPrice = yearFactor.*0.8.*GenPrice;

Create a new optimization problem named profitMaxExpandProb.

profitMaxExpandProb = optimproblem;

Set upper bounds on the installed capacity of each technology.

capacity.UpperBound = [100; 200; 150; 100; 300];

Define the profit to maximize, which includes the equivalent annual cost division factor.

weightedDispatchProfit = GenPrice./yearFactor - OperCost;
profitMaxExpandProb.Objective = sum(yearFactor.* ...
    sum( sum( weightedDispatchProfit .*loadDispatch, 3), 2) ) - ...
    dot(GenInvestCost,capacity)/divFactor;
profitMaxExpandProb.ObjectiveSense = "maximize";

Place the bounds on the capacity in the problem and solve the problem.

profitMaxExpandProb.Constraints.CapacityLimit = capLimit;
[solProfitExpand,totalProfit] = solve(profitMaxExpandProb)

Solving problem using linprog.

Optimal solution found.

solProfitExpand = struct with fields:
        capacity: [5×1 double]
    loadDispatch: [4×4×5 double]

totalProfit = 
7.1655e+07

View the optimal capacities of each technology.

solProfitExpand.capacity

ans = 5×1

   100
   200
   150
   100
     0

All the optimal capacities are at their upper limits except for nuclear, which is unused (0). View the dispatched power for each technology.

solProfitExpand.loadDispatch

ans = 
ans(:,:,1) =

  300.0000   60.0000  300.0000  300.0000
  120.0000  300.0000  450.0000  120.0000
  120.0000  300.0000  450.0000  120.0000
  300.0000  300.0000  450.0000  450.0000


ans(:,:,2) =

   1.0e+03 *

    0.6000    0.9000    0.6000    0.6000
    0.3600    0.9000    1.0200    0.6000
    0.3600    0.9000    1.0200    0.6000
    0.6000    0.6000    1.0200    1.0200


ans(:,:,3) =

         0  540.0000  540.0000         0
         0  630.0000  675.0000         0
         0  720.0000  900.0000  450.0000
         0  450.0000  675.0000  180.0000


ans(:,:,4) =

   600   600   600   600
     0     0     0     0
   600   600   600   600
     0     0     0     0


ans(:,:,5) =

     0     0     0     0
     0     0     0     0
     0     0     0     0
     0     0     0     0

Onshore and offshore wind (matrices 1 and 2) are used heavily, as is solar power (matrix 3) during periods when it is available. Gas (matrix 4) is used only in winter and summer, rows 1 and 3. Calculate the long-term average dispatched load.

avload = mean(solProfitExpand.loadDispatch,[1,2])

avload = 
avload(:,:,1) =

  277.5000


avload(:,:,2) =

  731.2500


avload(:,:,3) =

   360


avload(:,:,4) =

   300


avload(:,:,5) =

     0

In the problem units, each of the long-term average powers handled by onshore wind, solar, and gas is near 300. The average for offshore wind is over 700.

Conclusion

Deterministic problems of energy investment and power dispatch can be formulated and solved using linear programming. The problem-based approach eases the formulation of problems by making symbolic-style variables available.

You can generalize these models in several ways, such as adding ramping constraints that account for generators not being able to reach maximum capacity instantly, or having storage available for power with charging constraints and energy loss.