This article presents a fully-commented MATLAB implemantation that couples YALMIP modeling with CPLEX to reproduce the key findings of a recent SCI paper on bi-level, multi-timescale optimization of geographically-distributed energy hubs. The package is designed for rapid experimentation: every module is self-contained, vectorized, and ready for extension.

Modeling Highlights

The overall problem is decomposed into two decision layers that interact iteratively:

1. Local layer – each energy hub solves its own cost-minimization problem over a receding horizon while handling stochastic PV and load via a convex relaxation of the complementarity constraints.
2. Coordination layer – a central coordinator minimizes the aggregated operating cost of all hubs while enforcing tie-line capacity limits and preserving privacy (only incentive/response signals are exchanged).

To meet real-time requirements, a two-phase solution scheme is introduced: Phase-1 performs day-ahead scheduling, Phase-2 refines the schedule in a sliding-window fashion every 5–15 minutes.

Key Algorithmic Ingredients

• Convex relaxation: two mild sufficient conditions are derived to turn the non-convex complementarity constraints into second-order cone constraints, keeping the relaxation exact under practical operating ranges.
• Privacy-preserving coordination: a auxiliary problem principle (APP) variant is employed so that no sensitive parameters (cost curves, capacities) leave the local solvers.
• Iterative reduction: warm-starting Phase-2 with Phase-1 solutions typically cuts the iteration count by 60–70 % compared with a naïve fixed-point scheme.

Code Structure

The repository is organized into five functional blocks:

1. Price-Clearing Engine

function lambda = marketClear(bidStack, residualDemand)
% bidStack: monotonically decreasing vector of aggregated bids
% residualDemand: scalar, positive = buy, negative = sell
global priceMin priceMax deltaP

if residualDemand > bidStack(1)
    lambda = priceMin;
elseif residualDemand < bidStack(end)
    lambda = priceMax;
else
    idx = find(bidStack <= residualDemand, 1, 'first');
    lambda = priceMin + (idx-1)*deltaP;
end
end

2. Thermal Network Model

% Heat storage linearized dynamics
Aeq_T = [];  beq_T = [];
A1_Tsoc = zeros(H, nVar);  b1_Tsoc = zeros(H,1);
A2_Tsoc = zeros(H, nVar);  b2_Tsoc = zeros(H,1);

if hub.TES_cap > 0
    Aeq_T = zeros(1, nVar);
    beq_T = -(hub.TES_targetSOC - hub.TES_SOC * hub.TES_loss^H) * hub.TES_cap;
    for k = 1:H
        Aeq_T(1, idxDis(k)) = hub.TES_loss^(H-k)/hub.TES_eff;
        Aeq_T(1, idxChg(k)) = -hub.TES_loss^(H-k)*hub.TES_eff;
    end
end

3. Gas Network Parameters

% Gas import tariffs (converted to $/kWh)
gasTariff(1) = 0.334;   % peak rate
gasTariff(2) = 0.284;   % off-peak rate

% Flow limits (kWh)
gasCap(1) = 1e6;  gasCap(2) = 1e6;  gasCap(3) = 1e6;

4. Congession Management

% Tie-line power under LP vs MILP clearing
P_LP  = -sum(result(:,1:3),2);   % relaxed clearing
P_MILP= -sum(result(:,4:6),2);   % integer clearing

5. Demand-Response & Rolling Horizon

% Scenario flag: 0 = distributed, 2 = centralized
global runMode exportFlag
runMode   = 2;   % centralized
exportFlag= 1;   % allow energy export

initializeParams;          % load scenario data
offGrid = 0;              % 1 = island IES-1

% Day-ahead or intra-day switch
isDA = false;             % true = day-ahead, false = intra-day
singleWindowSolver;       % 5-min rolling update

Typical Outputs

Running the master script main.m produces:

• Inter-hub exchanges: 24-hour profile of tie-line flows and shadow prices.
• Local schedules: per-hub dispatch of CHP, boilers, PV curtailment, battery SOC, and heat storage SOC.

All figures are auto-generated and saved as vector PDFs for direct inclusion in LaTeX reports.

Extending the Package

The modular design allows quick swaps of:

• Network topologies (meshed vs radial)
• Uncertainty models (robust, stochastic, or distributionally robust)
• Market rules (pay-as-bid vs marginal pricing)

Each module is unit-tested; simply add new test cases under /tests.