Linear Programming for
Energy Systems
Theory, Geometry, and Implementation
ESE 4322 Energy Systems Optimization
What is Linear Programming?
A mathematical method to achieve the best outcome in a
mathematical model whose requirements are represented by linear
relationships.
Components:
1. Decision Variables (x)
2. Objective Function (Linear: fTx)
3. Constraints (Linear inequalities/equalities)
The 4 Assumptions of LP
1. Proportionality: Contribution to cost is proportional to value
(No economies of scale, x2).
2. Additivity: Total cost is sum of individual costs
(No interaction terms x*y).
3. Divisibility: Variables can be fractional
(Continuous, not Integer).
4. Certainty: Parameters (prices, limits) are known constants.
Geometry: Half-Spaces
Every linear inequality (Ax <= b) cuts the universe in half.
• Imagine a line on a 2D plot.
• One side is 'Valid' (Feasible).
• The other side is 'Invalid' (Infeasible).
• Adding more constraints chops away more space.
The Feasible Region
SThe intersection of all Half-Spaces.
It forms a Polytope (Polygon in n-
dimensions).
It contains EVERY valid solution to your
engineering problem.
If the region is empty? Your problem is
Infeasible.
Convexity
• LP Feasible Regions are always CONVEX.
• Definition: If you take any two points in the region, the line
connecting them is also inside the region.
• Implication: There are no 'local valleys' to get stuck in. A local
minimum is GUARANTEED to be the Global Minimum.
Visualizing the Objective
The Objective Function (min cTx) is a direction.
• Imagine 'Iso-Cost' lines perpendicular to the gradient.
• We 'slide' this line across the feasible region in the direction of
improvement.
• The last point we touch before leaving the region is the Optimum.
Types of LP Solutions
1. Infeasible: Constraints contradict (No region exists)
2. Unbounded: You can improve forever (Missing constraint).
3. Unique Optimal: Hits a single sharp corner
2. Multiple Optimal: The objective is parallel to an edge (Any point on
the edge is good)
The Fundamental Theorem of LP
• If an optimal solution exists, at least one optimal solution occurs
at a vertex (corner point) of the feasible region.
• Why? Because the center of the region is never better than the
edge in a linear world.
The Simplex Algorithm
Since the answer is at a corner, why search the inside?
Algorithm:
1. Start at ANY valid corner.
2. Look along the edges connected to this corner.
3. Does moving along an edge reduce cost?
4. If YES, move to next corner. Repeat.
5. If NO, you are at the Optimum. Stop.
Modern Solvers: Interior Point
Simplex walks the edges. Interior Point cuts
through the middle.
• For massive energy problems (1M+ variables),
Simplex can be slow.
• Interior Point (Barrier methods) are often
default in MATLAB.
• But the 'Vertex' concept remains the best way
to understand the solution.
Why Energy Fits LP
Most grid physics are naturally linear:
• Kirchhoff's Current Law (Sum flows = 0) → Linear Equality.
• Generator Limits (Min/Max MW) → Linear Bounds.
• Energy Storage (Bucket model) → Linear Dynamics.
• Cost ($/MWh) → Linear Objective.
Standard Form Transformation
MATLAB requires: Few useful tricks:
min fTx 1. Maximization? Minimize (-fTx).
2. Greater Than (>=)? Multiply by -1 to get
s.t.
(<=).
Ax <= b 3. Absolute Values? Split into two inequalities.
Aeq x = beq
The Engineering Workflow
1. Physics Model (Whiteboard): Pgrid(t) + Ppv(t) - Pbat(t) = Load
2. Mathematical Standard Form: Aeq * x = beq
3. Code (MATLAB/Python): x = linprog(f, A, b...)
Building the Matrix: Step 1 (Variables)
We must flatten time-series into a vector.
x = [ Pgrid(1..T) ; Ppv(1..T) ; Pbat(1..T) ; SOC(1..T) ]
• If T=4:
Indices 1-4: Grid
Indices 5-8: PV
Indices 9-12: Battery
Indices 13-16: SOC
Building Matrix: Step 2 (Power Balance)
• For every hour t, we have one equation:
Pgrid(t) + Ppv(t) - Pbat(t) = Load(t)
• This forms a 'Diagonal' structure in the Aeq matrix.
• It links variables of the *same* time step.
Building Matrix: Step 3 (Dynamics)
Battery SOC connects time steps:
SOC(t+1) - SOC(t) - Pbat(t) * dT = 0
This creates a 'Staircase' or 'Off-Diagonal' structure.
This is what makes the problem coupled over time.
Building Matrix: Step 4 (Limits)
Two ways to handle limits:
1. Simple Bounds (lb/ub vectors):
min <= Ppv <= max
2. Matrix Inequalities (A, b):
Used for complex limits like 'Ramping' (Pt - Pt-1 <= Limit).
Matrix Sparsity
• Real energy matrices are 99% Zeros.
• Variable at t=1 usually doesn't affect t=100 directly.
• Use MATLAB's 'sparse()' command to save memory.
• Spy plot visualizes the interaction structure.
MATLAB Implementation
• f = [price_vec; zeros; zeros; zeros]; % Objective
• Aeq = [...]; beq = [...]; % Physics
• lb = [...]; ub = [...]; % Limits
• [x, cost, exitflag] = linprog(f, [], [], Aeq, beq, lb, ub);
Practical Tip: Scaling
• Bad Scaling kills solvers.
• Don't mix Volts (10^0) with Grid Power (10^9 Watts).
• Normalize your units!
• Use MW (Megawatts) and PU (Per Unit) where possible.
Debugging Infeasibility
Exitflag = -2 (Infeasible).
What to do?
1. Check signs (did you swap <= and >= ?).
2. Check initial conditions (Is SOC_start within bounds?).
3. Slack Variables: Add a 'dummy generator' with high cost to see
where physics breaks.
Summary & Next Steps
1. LP is the bridge between Engineering Models and Optimal Decisions.
2. Convexity guarantees we find the best answer.
3. Simplex 'crawls the corners' of the feasible polygon.
Lab Session - We write this code from scratch.