Battery Energy Storage Optimization

🎯 Lab Objectives

• Understand battery scheduling as an optimization problem.
• Model the battery State of Charge.
• Include charging and discharging efficiencies.
• Respect SOC and power limits.
• Use electricity price as part of the objective function.
• Compare uncontrolled and optimized battery operation.
• Implement a simple optimization routine in MATLAB.

1. Battery Scheduling Problem

A battery can charge when electricity is cheap and discharge when electricity is expensive.

2. Time Horizon

Consider a six-hour scheduling horizon:

3. Battery Parameters

4. Battery Power Sign Convention

Define:

Therefore:

5. State of Charge

Battery energy is:

and:

6. SOC During Charging

If PB is negative, the battery is charging:

7. SOC During Discharging

If PB is positive:

The stored energy decreases slightly more than the delivered electrical energy because efficiency is less than one.

8. SOC Constraints

For this battery:

9. Grid Power

Without PV in this first example:

If the battery discharges, grid demand decreases.

If the battery charges, grid demand increases.

10. Electricity Cost

For one hour:

The total operating cost is:

11. Decision Vector

The optimization variables are battery powers over all hours:

One candidate schedule may be:

12. Why Do We Need Penalties?

A randomly generated schedule may violate SOC limits.

or:

Therefore define:

13. SOC Penalty

At every hour:

if SOC is too low, and:

if SOC is too high.

14. Final SOC Constraint

Without a final-energy condition, the optimizer may simply empty the battery near the end.

Therefore require:

and add:

15. MATLAB Data

price = [50 40 35 70 90 80];

loadMW = [6.0 5.5 5.0 6.5 7.0 6.5];

N = length(price);

dt = 1;

16. Battery Parameters in MATLAB

Emax = 4;

SOCmin = 0.20;
SOCmax = 0.90;
SOC0   = 0.50;

PchMax  = 2;
PdisMax = 2;

etaCh  = 0.95;
etaDis = 0.95;

17. Evaluate One Schedule

Start with:

Then simulate every hour sequentially.

E = SOC0*Emax;

for k = 1:N

    PB = schedule(k);

    if PB < 0

        E = E + etaCh*(-PB)*dt;

    else

        E = E - PB*dt/etaDis;

    end

    SOC = E/Emax;

end

18. Cost Calculation

Pgrid = loadMW(k) - PB;

cost = cost + price(k)*Pgrid*dt;

The cost is accumulated across the scheduling horizon.

19. Evolutionary Optimizer

We use a simple population-based optimizer.

20. Complete MATLAB Program

clear; clc; close all;

% ------------------------------------------------
% Battery Energy Storage Optimization
% Educational evolutionary scheduling example
% ------------------------------------------------

price = [50 40 35 70 90 80];

loadMW = [6.0 5.5 5.0 6.5 7.0 6.5];

N = length(price);

dt = 1;

% ------------------------------------------------
% Battery parameters
% ------------------------------------------------

Emax = 4;

SOCmin = 0.20;
SOCmax = 0.90;

SOC0 = 0.50;

PchMax  = 2;
PdisMax = 2;

etaCh  = 0.95;
etaDis = 0.95;

% ------------------------------------------------
% Optimization parameters
% ------------------------------------------------

popSize = 60;
maxGen  = 150;

Pc = 0.8;
Pm = 0.12;

sigmaMut = 0.4;

lambdaSOC   = 1e5;
lambdaFinal = 1e5;

% ------------------------------------------------
% Initial population
% Each row is one battery schedule
% ------------------------------------------------

population = ...
    -PchMax + ...
    rand(popSize,N)*(PchMax+PdisMax);

bestHist = zeros(maxGen,1);

% ------------------------------------------------
% Evolutionary loop
% ------------------------------------------------

for gen = 1:maxGen

    Jvec = zeros(popSize,1);

    for i = 1:popSize

        schedule = population(i,:);

        E = SOC0*Emax;

        totalCost = 0;
        penaltySOC = 0;

        for k = 1:N

            PB = schedule(k);

            % ------------------------------------
            % Battery energy update
            % ------------------------------------

            if PB < 0

                E = E + ...
                    etaCh*(-PB)*dt;

            else

                E = E - ...
                    PB*dt/etaDis;

            end

            SOC = E/Emax;

            % ------------------------------------
            % SOC penalties
            % ------------------------------------

            if SOC < SOCmin

                penaltySOC = penaltySOC + ...
                    lambdaSOC*(SOCmin-SOC)^2;

            elseif SOC > SOCmax

                penaltySOC = penaltySOC + ...
                    lambdaSOC*(SOC-SOCmax)^2;

            end

            % ------------------------------------
            % Grid power and electricity cost
            % ------------------------------------

            Pgrid = loadMW(k) - PB;

            totalCost = ...
                totalCost + ...
                price(k)*Pgrid*dt;

        end

        % ----------------------------------------
        % Final SOC penalty
        % ----------------------------------------

        SOCfinal = E/Emax;

        penaltyFinal = ...
            lambdaFinal*(SOCfinal-SOC0)^2;

        % ----------------------------------------
        % Total objective
        % ----------------------------------------

        Jvec(i) = ...
            totalCost + ...
            penaltySOC + ...
            penaltyFinal;

    end

    % ------------------------------------------------
    % Elitism
    % ------------------------------------------------

    [bestJ,bestIdx] = min(Jvec);

    elite = population(bestIdx,:);

    bestHist(gen) = bestJ;

    % ------------------------------------------------
    % New population
    % ------------------------------------------------

    newPopulation = zeros(size(population));

    newPopulation(1,:) = elite;

    row = 2;

    while row <= popSize

        % ----------------------------------------
        % Tournament selection - parent 1
        % ----------------------------------------

        a = randi(popSize);
        b = randi(popSize);

        if Jvec(a) < Jvec(b)
            parent1 = population(a,:);
        else
            parent1 = population(b,:);
        end

        % ----------------------------------------
        % Tournament selection - parent 2
        % ----------------------------------------

        a = randi(popSize);
        b = randi(popSize);

        if Jvec(a) < Jvec(b)
            parent2 = population(a,:);
        else
            parent2 = population(b,:);
        end

        % ----------------------------------------
        % Arithmetic crossover
        % ----------------------------------------

        if rand < Pc

            alpha = rand;

            child1 = ...
                alpha*parent1 + ...
                (1-alpha)*parent2;

            child2 = ...
                (1-alpha)*parent1 + ...
                alpha*parent2;

        else

            child1 = parent1;
            child2 = parent2;

        end

        % ----------------------------------------
        % Mutation
        % ----------------------------------------

        for k = 1:N

            if rand < Pm

                child1(k) = ...
                    child1(k) + ...
                    sigmaMut*randn;

            end

            if rand < Pm

                child2(k) = ...
                    child2(k) + ...
                    sigmaMut*randn;

            end

        end

        % ----------------------------------------
        % Battery power limits
        % ----------------------------------------

        child1 = ...
            min(max(child1,-PchMax),PdisMax);

        child2 = ...
            min(max(child2,-PchMax),PdisMax);

        % ----------------------------------------
        % Store offspring
        % ----------------------------------------

        newPopulation(row,:) = child1;

        if row+1 <= popSize
            newPopulation(row+1,:) = child2;
        end

        row = row + 2;

    end

    population = newPopulation;

end

% ------------------------------------------------
% Final evaluation
% ------------------------------------------------

Jvec = zeros(popSize,1);

for i = 1:popSize

    schedule = population(i,:);

    E = SOC0*Emax;

    totalCost = 0;
    penaltySOC = 0;

    for k = 1:N

        PB = schedule(k);

        if PB < 0

            E = E + ...
                etaCh*(-PB)*dt;

        else

            E = E - ...
                PB*dt/etaDis;

        end

        SOC = E/Emax;

        if SOC < SOCmin

            penaltySOC = penaltySOC + ...
                lambdaSOC*(SOCmin-SOC)^2;

        elseif SOC > SOCmax

            penaltySOC = penaltySOC + ...
                lambdaSOC*(SOC-SOCmax)^2;

        end

        Pgrid = loadMW(k) - PB;

        totalCost = ...
            totalCost + ...
            price(k)*Pgrid*dt;

    end

    SOCfinal = E/Emax;

    penaltyFinal = ...
        lambdaFinal*(SOCfinal-SOC0)^2;

    Jvec(i) = ...
        totalCost + ...
        penaltySOC + ...
        penaltyFinal;

end

[bestJ,bestIdx] = min(Jvec);

bestSchedule = population(bestIdx,:);

% ------------------------------------------------
% Simulate best schedule
% ------------------------------------------------

SOCbest = zeros(1,N+1);
PgridBest = zeros(1,N);

SOCbest(1) = SOC0;

E = SOC0*Emax;

bestCost = 0;

for k = 1:N

    PB = bestSchedule(k);

    if PB < 0

        E = E + ...
            etaCh*(-PB)*dt;

    else

        E = E - ...
            PB*dt/etaDis;

    end

    SOCbest(k+1) = E/Emax;

    PgridBest(k) = ...
        loadMW(k) - PB;

    bestCost = ...
        bestCost + ...
        price(k)*PgridBest(k)*dt;

end

% ------------------------------------------------
% Cost without battery
% ------------------------------------------------

baseCost = ...
    sum(price .* loadMW * dt);

saving = baseCost - bestCost;

% ------------------------------------------------
% Display results
% ------------------------------------------------

fprintf('Battery Optimization Result\n\n')

fprintf('Base cost without battery = %.2f\n', ...
        baseCost)

fprintf('Optimized cost = %.2f\n', ...
        bestCost)

fprintf('Cost saving = %.2f\n\n', ...
        saving)

fprintf('Battery schedule [MW]:\n')
disp(bestSchedule)

fprintf('SOC trajectory:\n')
disp(SOCbest)

% ------------------------------------------------
% Plot convergence
% ------------------------------------------------

figure

plot(1:maxGen,bestHist,'LineWidth',1.6)

grid on

xlabel('Generation')
ylabel('Best Objective')

title('Battery Scheduling Optimization Convergence')

% ------------------------------------------------
% Plot battery power
% ------------------------------------------------

figure

bar(1:N,bestSchedule)

grid on

xlabel('Hour')
ylabel('Battery Power [MW]')

title('Optimal Battery Schedule')

% ------------------------------------------------
% Plot SOC
% ------------------------------------------------

figure

plot(0:N,SOCbest,'o-','LineWidth',1.6)

hold on

yline(SOCmin,'--')
yline(SOCmax,'--')

grid on

xlabel('Hour')
ylabel('State of Charge')

title('Battery SOC Trajectory')

% ------------------------------------------------
% Plot electricity price
% ------------------------------------------------

figure

stairs(1:N,price,'LineWidth',1.6)

grid on

xlabel('Hour')
ylabel('Electricity Price [EUR/MWh]')

title('Electricity Price Profile')

21. Cost Without the Battery

The reference operating cost is:

MATLAB:

baseCost = sum(price .* loadMW * dt);

22. Economic Saving

A positive value means the optimized battery schedule reduced electricity purchasing cost.

23. Interpreting the Optimal Schedule

24. Experiment 1 — Change Initial SOC

SOC0 = 0.30;
SOC0 = 0.50;
SOC0 = 0.80;

Compare optimal schedules and savings.

25. Experiment 2 — Change Battery Capacity

Emax = 2;
Emax = 4;
Emax = 8;

26. Experiment 3 — Change Efficiency

etaCh  = 1.00;
etaDis = 1.00;

Then compare with:

etaCh  = 0.90;
etaDis = 0.90;

27. Experiment 4 — Flat Electricity Price

Try:

price = [60 60 60 60 60 60];

Ask whether battery arbitrage still provides an economic benefit.

28. Add PV Generation

A natural extension is:

The battery can then store excess PV energy and use it later.

29. Battery Degradation Cost

Real optimization should also consider battery aging.

A simple educational approximation is:

Then:

30. From Scheduling to Smart Energy Management

A more advanced energy-management problem may optimize:

and minimize:

31. Lab Tasks

1. Enter the electricity-price and load profiles.
2. Implement charging and discharging equations.
3. Calculate SOC hour by hour.
4. Add SOC penalties.
5. Add the final-SOC condition.
6. Run the evolutionary optimizer.
7. Plot the optimal battery schedule.
8. Plot the SOC trajectory.
9. Calculate the cost without a battery.
10. Calculate the optimized cost and savings.
11. Change the initial SOC.
12. Change battery capacity and efficiency.
13. Test a flat electricity-price profile.
14. Explain how PV could be added to the model.

📝 Review Questions

1. What does positive battery power mean in this laboratory?
2. How does charging efficiency affect SOC?
3. Why is SOC calculated sequentially?
4. Why do we need SOC limits?
5. Why is a final-SOC constraint useful?
6. Why does the battery tend to charge during low-price hours?
7. How does low efficiency affect economic savings?
8. Why may a flat electricity price reduce the value of storage?
9. Why should degradation cost be included in more realistic studies?
10. How would PV generation modify the grid-power equation?