Computational burdoen of bilevel problem
20 views
Skip to first unread message
weiqi meng
Jun 7, 2025, 1:32:23 PMJun 7
to YALMIP
Hi,
I just use kkt() to solve bilevel problem. but the computationa time is really heavy. Th inner level problem is just an LP problem. Attached is the code. Does tthere better way to improve the computatinal efficiency?
% function result_SP = test_123_kktsystem_SP(data_MP)
%% Reduce redundant defination of variables in KKT system.
Sbase = 100; % 100 kVA
Vbase = 4.16;
%% Load data from MP
% p_pv = data_MP.p_pv; % PV scenarios (n_der × 24 × Iterations)
% [n_der, n_time, Iterations, priceCT] = size(p_pv); % Number of iterations (scenarios) from subproblem
%
% p_DA = data_MP.p_DA; % 24*10 dimension constants
%% check use
p_DA = [20 3.62622428250000 3.62622428250000 3.62622428250000 3.62622428250000 3.62622428250000 3.62622428250000 3.62622428250000 3.58235078250000 3.62622428250000
9.17704460500002 -10.8229553950000 -20 -20 -10.8229553950000 -10.8229553950000 -20 -10.8229553950000 -20 -20
9.07292388499996 -10.9270761150000 -20 -20 -10.9270761150000 -10.9270761150000 -20 -10.9270761150000 -20 -20
-4.10863844303328 -10.9645233350000 -10.9645233350000 -10.9645233350000 -10.9645233350000 -10.9645233350000 -10.9645233350000 -10.9645233350000 -10.9645233350000 -10.9645233350000
-3.12863563215369 -3.12863563215369 -3.12863563215369 -3.63517383473750 -3.12863563215369 -3.12863563215369 -17.4422793233033 -3.12863563215369 -3.12863563215369 -17.4009640878462
8.81262780500001 -11.1873721950000 -11.1873721950000 -20 -20 -11.1873721950000 -20 -20 -20 -20
20 8.76061763249996 7.07725962249993 -4.74659769299998 -4.74659769299998 8.76061763249996 -4.74659769299998 -4.74659769299998 -11.5328423675000 -4.74659769299998
20 8.70860746000002 -11.2913925400000 -11.2913925400000 -11.2913925400000 8.70860746000002 -11.2913925400000 -11.2913925400000 -11.2913925400000 -11.2913925400000
20 20 -3.26519376000007 -3.26519376000007 -3.26519376000007 20 -4.15417176000000 -3.26519376000007 -3.26519376000007 -4.15417176000000
20 20 10.6763660200000 10.6763660200000 20 20 10.6763660200000 10.6763660200000 10.6763660200000 10.6763660200000
20 20 20 20 20 20 19.4420309725000 20 18.5083074996412 20
20 11.2394323755000 20 11.2394323755000 20 11.2394323755000 11.2394323755000 11.2394323755000 6.98431270250002 11.2394323755000
20 20 20 20 20 20 20 20 10.6066715200000 10.5646280200000
20 20 20 3.85887824000001 20 20 3.85887824000001 3.85887824000001 -16.2438807600000 -16.2438807600000
20 20 20 20 20 20 -15.6913114040997 20 20 20
20 20 20 20 20 20 6.84593682000000 20 20 11.8757109557313
20 20 20 18.9032544925000 20 20 18.9032544925000 20 20 -20
20 20 20 -20 -3.13707920250000 20 -11.1873721950000 20 -11.1873721950000 -11.1873721950000
20 8.70860746000002 8.70860746000002 8.63543596000002 8.63543596000002 8.70860746000002 -20 8.70860746000002 8.63543596000002 8.63543596000002
14.0718942775000 -11.4996339800000 8.50036601999997 -11.4996339800000 -11.4996339800000 8.50036601999997 -11.4996339800000 8.50036601999997 -11.4996339800000 -11.4996339800000
8.60448673999996 8.60448673999996 8.60448673999996 -11.3955132600000 -11.3955132600000 8.60448673999996 -11.3955132600000 8.60448673999996 -11.3955132600000 -11.3955132600000
8.81262780500001 8.81262780500001 8.81262780500001 8.81262780500001 8.81262780500001 8.81262780500001 -17.2406447711772 8.81262780500001 -17.2406447711772 -17.2406447711772
-3.59378803862881 -3.59378803862881 9.02088871500001 -10.9791112850000 -3.59378803862881 -3.59378803862881 -10.9791112850000 -3.59378803862881 -10.9791112850000 -10.9791112850000
5.71337039915513 5.71337039915513 5.71337039915513 2.95250040452910 2.95250040452910 5.71337039915513 -10.8749905650000 2.95250040452910 -10.8749905650000 -10.8749905650000];
% p_DA = p_DA';
DER_BusName = {'S1a';'S1b';'S1c';'S29a';'S29b';'S29c';'S35a';'S35b';'S35c';
'S44a';'S44b';'S44c';'S51a';'S51b';'S51c';'S86a';'S86b';'S86c';
'S108a';'S108b';'S108c'};
n_der = length(DER_BusName); % n_der = 21
% PV initial uncertainty distribution (1×24)
p_pv_base = 4*5 * [0 0 0 0 0 0 0 0 0.6 1.76 1.76 1.5 1.78 1.53 1.48 1.46 0.81 0 0 0 0 0 0 0];
p_pv_base = p_pv_base / Sbase;
% Initialize data_SP.p_pv with n_der × 24 × Iterations × priceCT
data_MP.p_pv = repmat(reshape(p_pv_base, [1, 24, 1, 1]), n_der, 1, 1, 10); % n_der × 24 × 1 × 10
p_pv = data_MP.p_pv; % PV scenarios (n_der × 24 × Iterations)
[n_der, n_time, priceCT, Iterations] = size(p_pv); % Number of iterations (scenarios) from subproblem
%% Load price scenario data (day-ahead and real-time)
load reduced_da_scenarios.mat
load reduced_rt_scenarios.mat
load scenario_probabilities.mat
lambda_da = Sbase * reduced_da_scenarios' /1000; % 10*24价格矩阵 ($/MWh)
rho_da = finalProbs; % 10x1 probability vector
lambda_rt = Sbase * reduced_rt_scenarios' /1000; % 10*24价格矩阵 ($/MWh)
rho_rt = finalProbs; % 10x1 probability vector
priceCT = size(rho_rt, 1); % Number of scenarios (10)
% priceCT = 10;
% load DLMP_data
% load DLMP_data_RT
% a_data=[a_data(9:24);a_data(1:8)];b_data=[b_data(9:24);b_data(1:8)];price_RT=[price_RT(:,9:24),price_RT(:,1:8)]';
% a_data = repmat(a_data,1,priceCT);
% b_data = repmat(b_data,1,priceCT);
% lambda_da = b_data;
% k_lambda_da = a_data;
% lambda_rt = price_RT;
% rho_rt = PDF;
% rho_da = PDF;
% priceCT = 10;
% load DLMP_data
% load DLMP_data_RT
% a_data=[a_data(9:24);a_data(1:8)];b_data=[b_data(9:24);b_data(1:8)];price_RT=[price_RT(:,9:24),price_RT(:,1:8)]';
% a_data = repmat(a_data,1,priceCT); % 24*10
% b_data = repmat(b_data,1,priceCT);
% lambda_rt = b_data;
% k_lambda_rt = a_data;
% lambda_da = price_RT;
% rho_rt = PDF;
% rho_da = PDF;
% Load network data
bus_data = load('load_name_bus_phase.mat');
load_name_bus_phase = bus_data.load_name_bus_phase;
ldf_data = load('123_LDF.mat');
A = ldf_data.A;
A0 = ldf_data.A0;
RD = ldf_data.RD;
XD = ldf_data.XD;
load_data = load('123_load.mat');
load_Name = load_data.load_Name;
load_kW = load_data.load_kW;
load_kVar = load_data.load_kVar;
% Number of buses and lines
nb = size(A,1);
nl = nb;
% DER bus indices
DER_BusName = {'S1a';'S1b';'S1c';'S29a';'S29b';'S29c';'S35a';'S35b';'S35c';
'S44a';'S44b';'S44c';'S51a';'S51b';'S51c';'S86a';'S86b';'S86c';
'S108a';'S108b';'S108c'};
DER_bus = zeros(length(DER_BusName),1);
for i = 1:length(DER_BusName)
idx = find(strcmp(load_name_bus_phase(:,1), DER_BusName{i}), 1);
DER_bus(i) = idx;
end
n_der = length(DER_bus);
% Load data (pu, nb x 24)
p_load = zeros(nb, 24);
q_load = zeros(nb, 24);
% 24-hour load distribution factors
probabilities = [
0.026503717, 0.028398378, 0.032184586, 0.035970794, 0.037862073, 0.041649898, ...
0.043541177, 0.045432456, 0.049218664, 0.053004872, 0.054896151, 0.056787429, ...
0.053004872, 0.049218664, 0.045432456, 0.039753352, 0.037862073, 0.041649898, ...
0.045432456, 0.053004872, 0.049218664, 0.041649898, 0.034076774, 0.030290566];
for i = 1:size(load_Name, 1)
idx = find(strcmp(load_name_bus_phase(:,1), load_Name(i)), 1);
for t = 1:24
p_load(idx, t) = load_kW(i) * probabilities(t) / Sbase;
q_load(idx, t) = load_kVar(i) * probabilities(t) / Sbase;
end
end
% PV data (24-hour profile, kW)
p_pv_0 = 12 * 5 * [0 0 0 0 0 0 0 0 0.6 1.76 1.76 1.5 1.78 1.53 1.48 1.46 0.81 0 0 0 0 0 0 0];
% PV uncertainty bounds
delta = 0.5 * max(p_pv_0);
p_pv_real = repmat(p_pv_0 / Sbase, n_der, 1); % n_der × 24
p_pv_min = zeros(n_der, 24);
p_pv_max = zeros(n_der, 24);
non_zero = p_pv_0 > 0;
p_pv_min(:, non_zero) = repmat(max(0, p_pv_0(non_zero) - delta) / Sbase, n_der, 1);
p_pv_max(:, non_zero) = repmat((p_pv_0(non_zero) + delta) / Sbase, n_der, 1);
% Fload data (n_der x 24)
% Original 7 x 24 Fload_0 data for specific DER buses
Fload_0_orig = 0.05 * [28.403 27.745 25.595 28.858 27.907 38.154 35.231 41.2 33.813 37.836 44.419 43.83 41.183 40.113 38.614 44.702 52.482 53.223 44.659 42.954 42.16 33.82 27.996 25.726;
26.708 24.858 25.325 26.854 27.985 39.732 38.573 38.658 32.45 30.887 30.881 37.759 41.552 34.551 40.839 50.903 50.605 46.766 51.589 49.475 41.592 32.965 29.287 24.987;
26.12 23.628 25.965 23.224 24.335 28.837 36.887 34.475 34.442 35.13 29.114 33.104 32.029 36.393 30.638 38.503 48.347 46.176 46.177 44.766 40.912 35.985 27.62 26.049;
25.179 23.709 23.727 29.5 28.324 31.996 37.915 40.539 33.414 34.971 36.674 39.532 33.595 37.064 35.727 47.89 52.252 50.777 53.641 55.291 43.828 42.018 33.738 29.795;
33.249 26.929 27.074 26.424 30.609 39.057 34.964 36.167 48.314 39.416 42.604 41.926 45.965 40.968 41.437 45.023 47.052 43.705 36.877 39.688 38.824 37.688 32.84 28.788;
31.949 27.441 26.953 26.404 30.105 33.339 35.819 40.983 37.922 42.707 46.129 48.732 48.593 46.93 49.832 57.155 55.125 52.432 54.479 45.009 36.146 34.71 35.04 30.189;
27.27 26.902 26.891 24.765 27.193 30.024 35.482 45.443 49.897 53.75 50.579 60.073 64.758 56.226 68.008 61.09 59.055 57.432 52.781 56.064 45.286 39.492 33.21 29.44];
% Map to n_der × 24
Fload_0 = kron(Fload_0_orig, ones(3,1)); % n_der × 24
Fload_0_b = Fload_0 / Sbase;
Fload_adj = Fload_0 - 0.2;
Fload_max_kW = Fload_adj; % n_der × 24
Fload_min_kW = zeros(n_der, 24);
% ESS parameters
N_es = 10; % Number of Energy
E_ess_0 = N_es * 14.5; % 14.5 kWh = 0.145 pu-h
P_ess_max_0 = N_es * 11.3; % 11.3 kW = 0.113 pu
eta = 0.95; % Round-trip efficiency (sqrt for charge/discharge)
M = N_es * 11.3 / Sbase; % Big-M constant (sufficiently large)
% Cost coefficients ($/pu)
costVec = [1; 1.5; 1.2; 2; 1.8; 1.4; 2.2];
P_Fload_cost = kron(costVec, ones(3, 24)); % n_der × 24
PV_cost = 0.01 * ones(n_der, 1); % n_der × 1
ESS_cost = 0.1 * ones(n_der, 1); % n_der × 1
% Parameters
V_min = 0.9; V_max = 1.10;
v_min = V_min^2; v_max = V_max^2;
v0 = [1.00^2; 1.00^2; 1.00^2];
v0_term = -repmat(inv(A.') * A0 * v0, 1, 24);
p_Fload_max = Fload_max_kW / Sbase; % n_der × 24
p_Fload_min = Fload_min_kW / Sbase; % n_der × 24
P_ess_max = P_ess_max_0 * ones(n_der, 24) / Sbase; % n_der × 1
E_ess = E_ess_0 * ones(n_der, 24) / Sbase; % n_der × 1
Lambda = 0.5; % Uncertainty budget
% Initialize result arrays for all scenarios
p_Fload_all = zeros(n_der, 24, priceCT);
p_pv_all = zeros(n_der, 24, priceCT);
p_ess_ch_all = zeros(n_der, 24, priceCT);
p_ess_dis_all = zeros(n_der, 24, priceCT);
soc_all = zeros(n_der, 24, priceCT);
p_all = zeros(nb, 24, priceCT);
p_RT = zeros(24, priceCT);
q_all = zeros(nb, 24, priceCT);
p_slack_all = zeros(1, 24, priceCT);
q_slack_all = zeros(1, 24, priceCT);
obj_value_all = zeros(priceCT, 1);
result_info = cell(priceCT, 1);
slack_pos_all = zeros(24,priceCT); % Positive slack
slack_neg_all = zeros(24,priceCT); % Negative slack
% Parallel loop over scenarios
parfor s = 1:priceCT
% Define variables for scenario s
p_Fload = sdpvar(n_der, 24, 'full'); % DER active power
p_pv = sdpvar(n_der, 24, 'full'); % PV active power
p_ess_ch = sdpvar(n_der, 24, 'full'); % ESS charge power (positive)
p_ess_dis = sdpvar(n_der, 24, 'full'); % ESS discharge power (positive)
b_ess = sdpvar(n_der, 24, 'full'); % Binary variable for charge/discharge
soc = sdpvar(n_der, 24, 'full'); % ESS state of charge
p = sdpvar(nb, 24, 'full'); % Net active power
q = sdpvar(nb, 24, 'full'); % Net reactive power
v = sdpvar(nb, 24, 'full'); % Voltage squared
p_RT_s = sdpvar(24, 1, 'full'); % Real-time power deviation
slack_pos = sdpvar(24, 1, 'full'); % Positive slack
slack_neg = sdpvar(24, 1, 'full'); % Negative slack
% Lower-level constraints
constraints_lower = [];
% Define non-DER buses
all_buses = 1:nb;
non_DER_bus = setdiff(all_buses, DER_bus);
% Matrix-based constraints for all time steps (t = 1 to 24)
constraints_lower = [constraints_lower;
p(non_DER_bus,:) == p_load(non_DER_bus,:); % Non-DER buses: p = p_load
p(DER_bus,:) == p_load(DER_bus,:) + (Fload_0_b - p_Fload) - p_pv - p_ess_dis + p_ess_ch; % DER buses
q == q_load; % Reactive power balance
v == v0_term - 2 * RD * p - 2 * XD * q; % Voltage equation
p_Fload >= p_Fload_min; % Flexible load lower bound
p_Fload <= p_Fload_max; % Flexible load upper bound
p_pv >= p_pv_min; % PV lower bound
p_pv <= p_pv_max; % PV upper bound
v >= v_min; % Voltage lower bound
v <= v_max; % Voltage upper bound
p_ess_ch >= 0; % ESS charge power non-negativity
p_ess_ch <= P_ess_max; % ESS charge power upper bound
p_ess_dis >= 0; % ESS discharge power non-negativity
p_ess_dis <= P_ess_max; % ESS discharge power upper bound
p_ess_ch <= M * b_ess; % Big-M: charge if b_ess=1
p_ess_dis <= M * (1 - b_ess); % Big-M: discharge if b_ess=0
b_ess >= 0; b_ess <= 1;
soc >= 0; % SOC non-negativity
soc <= E_ess]; % SOC upper bound
% SOC dynamics (requires loop due to temporal coupling)
for t = 1:24
if t == 1
constraints_lower = [constraints_lower;
soc(:,t) == 0.5 * E_ess(:,t) + p_ess_ch(:,t) * eta - p_ess_dis(:,t) / eta];
else
constraints_lower = [constraints_lower;
soc(:,t) == soc(:,t-1) + p_ess_ch(:,t) * eta - p_ess_dis(:,t) / eta];
end
end
% Final SOC constraint
constraints_lower = [constraints_lower; soc(:,24) >= 0.5 * E_ess(:,t)];
% No charge/discharge constraint
neg_price_idx = find(lambda_rt(s,:) < 0);
if ~isempty(neg_price_idx)
constraints_lower = [constraints_lower; p_ess_ch(:,neg_price_idx) <= 0.2 * P_ess_max(:,neg_price_idx)];
end
%% Add the RT power constraint here
% Add the RT power constraint using p_RT_s (only definition)
constraints_lower = [constraints_lower;
p_RT_s == (-sum(p, 1)' - p_DA(:,s)); % Define p_RT_s
];
% Lower-level objective (maximize profit for scenario s, with uniform time weighting)
obj_lower = sum(lambda_rt(:,s) .* p_RT_s) - ...
sum(sum(P_Fload_cost .* p_Fload)) - ... % DER node power cost
ESS_cost' * sum(p_ess_dis, 2) - ... % ESS discharge cost
1e3 * sum(max(p_RT_s - 10, 0) + max(-p_RT_s - 10, 0)); % Penalty for p_RT_s > 10 or p_RT_s < -10
% obj_lower = - ( sum(lambda_rt(s,:) .* sum(p, 1)) + ...% Slack bus power cost for scenario s
% sum(sum(P_Fload_cost .* p_Fload)) + ... % DER node power cost
% PV_cost' * sum(p_pv, 2) + ... % PV node power cost
% ESS_cost' * sum(p_ess_dis, 2) ); % ESS discharge cost
% Upper-level objective: Minimize cost for scenario s
obj_upper = - ( sum(lambda_rt(:,s) .* p_RT_s ) + ...
sum(lambda_da(:,s) .* p_DA(:,s)) ) + ...% Slack bus power cost for scenario s
sum(sum(P_Fload_cost .* p_Fload)) + ... % DER node power cost
PV_cost' * sum(p_pv, 2) + ... % PV node power cost
ESS_cost' * sum(p_ess_dis, 2); % ESS discharge cost
% Upper-level constraints (PV uncertainty)
constraints_upper = [];
for j = 1:n_der
constraints_upper = [constraints_upper;
p_pv(j,:) >= p_pv_min(j,:);
p_pv(j,:) <= p_pv_max(j,:)];
% Uncertainty budget
sum_deviation = sum(p_pv_max(j,:) - p_pv(j,:));
sum_range = sum(p_pv_max(j,:) - p_pv_min(j,:));
if sum_range > 0
constraints_upper = [constraints_upper;
sum_deviation / sum_range <= Lambda];
end
end
% KKT system for lower-level problem
ops_kkt = sdpsettings('kkt.dualbound',0);
% ops_kkt = sdpsettings('kkt.dualbound',0, 'gurobi.FeasibilityTol',1e-9, 'gurobi.OptimalityTol',1e-9);
[kkt_constr, ~] = kkt(constraints_lower, -obj_lower, p_pv(:), ops_kkt);
% Single-level constraints
constraints = [constraints_upper; kkt_constr];
% Solve single-level problem
ops = sdpsettings('solver','gurobi','verbose',2, 'gurobi.Heuristics', 0.5, 'gurobi.TuneTimeLimit', 0, ...
'gurobi.NodefileStart', 2.0, 'gurobi.Threads', 4, 'gurobi.MIPGap', 0.001, ...
'gurobi.Heuristics', 0.05, ...
'gurobi.TimeLimit', 600, ...
'gurobi.MIPFocus', 2, 'gurobi.Presolve', 2,'showprogress',1); % MILP solver settings
% 'gurobi.MIPGap', 0.1, ... 'gurobi.TimeLimit', 600, ...
% 'gurobi.Threads', 10, ... 'gurobi.PreSparsify', 1, ... 'gurobi.PreSparsify', 1, ...
% ops = sdpsettings('solver','gurobi','verbose',3, 'gurobi.TuneTimeLimit',0, ...
% 'gurobi.MIPFocus',1, 'gurobi.Presolve',1, 'gurobi.FeasibilityTol',1e-9, ...
% 'gurobi.OptimalityTol',1e-9, 'showprogress',1); % MILP solver settings
% Optimize for scenario s
result = optimize(constraints, -obj_upper, ops);
% Store results for scenario s
p_Fload_all(:,:,s) = double(p_Fload); % n_der × 24
p_pv_all(:,:,s) = double(p_pv); % n_der × 24
p_ess_ch_all(:,:,s) = double(p_ess_ch); % n_der × 24
p_ess_dis_all(:,:,s) = double(p_ess_dis); % n_der × 24
soc_all(:,:,s) = double(soc); % n_der × 24
p_all(:,:,s) = double(p); % nb × 24
q_all(:,:,s) = double(q); % nb × 24
p_slack_all(1,:,s) = double(sum(p,1));
q_slack_all(1,:,s) = double(sum(q,1));
obj_value_all(s) = double(-obj_upper);
result_info{s} = result;
p_RT(:,s)=double(p_RT_s);
slack_pos_all=double(slack_pos);
slack_neg_all=double(slack_neg);
end
% Compute expected cost weighted by scenario probabilities
expected_profit = sum(rho_rt .* obj_value_all);
% Store results in output structure
result_SP.p_Fload = p_Fload_all;
result_SP.p_DA = p_DA;
result_SP.p_RT = p_RT;
% result_SP.p_pv = p_pv_all; % n_der × 24 × priceCT
result_SP.p_pv = reshape(p_pv_all, [n_der, 24, 1, priceCT]); % n_der × 24 × 1 × priceCT
result_SP.p_ess_ch = p_ess_ch_all;
result_SP.p_ess_dis = p_ess_dis_all;
result_SP.soc = soc_all;
result_SP.p = p_all;
result_SP.q = q_all;
result_SP.p_slack = p_slack_all;
result_SP.q_slack = q_slack_all;
result_SP.obj = expected_profit;
result_SP.result_info = result_info;
result_SP.obj_value_all = obj_value_all;
result_SP.p_load = sum(sum(p_load));
result_SP.p_load_sum = p_load(DER_bus,:);
result_SP.p_Fload_0 = Fload_0_b;
result_SP.p_Fload_0_sum = sum(sum(Fload_0_b));
result_SP.p_pv_1_sum= sum(sum(p_pv_all(:,:,1)));
result_SP.p_Fload_1_sum= sum(sum(Fload_0_b)) - sum(sum(p_Fload_all(:,:,1)));
result_SP.p_ess_ch_1_sum = sum(sum(p_ess_ch_all(:,:,1)));
result_SP.p_ess_dis_1_sum = sum(sum(p_ess_dis_all(:,:,1)));
result_SP.p_RT_sum = sum(sum(p_RT(:,:,1)));
result_SP.p_DA_sum = sum(sum(p_DA));
% Display results
disp('Bilevel problem solved for all scenarios!');
disp('Expected cost ($):');
disp(expected_profit); % Display as positive cost
for s = 1:priceCT
fprintf('Scenario %d (probability %.4f):\n', s, rho_rt(s));
if result_info{s}.problem == 0
fprintf(' Optimal cost ($): %.2f\n', -obj_value_all(s));
fprintf(' Optimal DER power output (pu, first 5 DER buses, all hours):\n');
disp(p_Fload_all(1:min(5,n_der),:,s));
fprintf(' Optimal PV power output (pu, first 5 DER buses, all hours):\n');
disp(p_pv_all(1:min(5,n_der),:,s));
fprintf(' Optimal ESS charge power (pu, first 5 DER buses, all hours):\n');
disp(p_ess_ch_all(1:min(5,n_der),:,s));
fprintf(' Optimal ESS discharge power (pu, first 5 DER buses, all hours):\n');
disp(p_ess_dis_all(1:min(5,n_der),:,s));
fprintf(' Optimal ESS SOC (pu, first 5 DER buses, all hours):\n');
disp(soc_all(1:min(5,n_der),:,s));
fprintf(' Slack bus active power (pu, all hours):\n');
disp(p_slack_all(1,:,s));
fprintf(' Slack bus reactive power (pu, all hours):\n');
disp(q_slack_all(1,:,s));
else
fprintf(' Problem not solved. Check solver output for details.\n');
disp(result_info{s}.info);
end
end
% end
Johan Löfberg
Jun 8, 2025, 2:42:15 AMJun 8
to YALMIP
No, you would have to develop your own dedicated solver. The MILP is simply too large and too weak.
Reply all
Reply to author
Forward
0 new messages