MATLAB-YALMIP optimal trendy program cannot converge
I use MATLAB-YALMIP to write the optimal trend program based on the IEEEE14 system, but I can’t converge. Please help to see if the program is wrong?
code:
% Optimal Power Flow Example
% Based on the IEEE 14-bus standard case (i.e., case14 in Matpower)
% Implemented in MATLAB with YALMIP
clear; clc;
% Load power system data
mpc = loadcase(‘case14’); % You can also load other data files, such as case9, case30, etc.
% Define variables and parameters
nb = size(mpc.bus, 1); % Number of buses (nodes)
nl = size(mpc.branch, 1); % Number of branches (lines)
ng = size(mpc.gen, 1); % Number of generators
nr = nb – ng; % Number of load buses
% Extract node data
bus_id = mpc.bus(:, 1); % Bus IDs
bus_type = mpc.bus(:, 2); % Bus types (1-PQ, 2-PV, 3-Slack)
Pd = mpc.bus(:, 3); % Active power demand (MW)
Qd = mpc.bus(:, 4); % Reactive power demand (MVAr)
Vmin = mpc.bus(:, 13); % Minimum bus voltage (p.u.)
Vmax = mpc.bus(:, 12); % Maximum bus voltage (p.u.)
% Extract generator data
gen_id = mpc.gen(:, 1); % Bus IDs where generators are connected
Pgmin = mpc.gen(:, 10); % Minimum generator active power output (MW)
Pgmax = mpc.gen(:, 9); % Maximum generator active power output (MW)
Qgmin = mpc.gen(:, 5); % Minimum generator reactive power output (MVAr)
Qgmax = mpc.gen(:, 4); % Maximum generator reactive power output (MVAr)
a2 = mpc.gencost(:, 5); % Generator cost curve coefficient for quadratic term
a1 = mpc.gencost(:, 6); % Generator cost curve coefficient for linear term
a0 = mpc.gencost(:, 7); % Generator cost curve constant term
% Extract branch data
fbus = mpc.branch(:, 1); % "From" bus (starting node) of each branch
tbus = mpc.branch(:, 2); % "To" bus (ending node) of each branch
Plmax = ones(nl, 1) * 4000; % Maximum branch active power flow limit (MW)
% Plmax = mpc.branch(:, 6) .* mpc.branch(:, 3) * 1e-3; % Maximum branch active power flow limit (MW)
% Calculate the bus admittance matrix
[Ybus, ~, ~] = makeYbus(mpc); % Use a Matpower function to calculate the bus admittance matrix
Gbus = real(Ybus); % Real part of bus admittance matrix (conductance)
Bbus = imag(Ybus); % Imaginary part of bus admittance matrix (susceptance)
% Define YALMIP variables
Pg = sdpvar(ng, 1, ‘full’); % Generator active power output variables
Qg = sdpvar(ng, 1, ‘full’); % Generator reactive power output variables
V = sdpvar(nb, 1, ‘full’); % Bus voltage magnitude variables
theta = sdpvar(nb, 1, ‘full’); % Bus voltage angle variables
Pl = sdpvar(nl, 1, ‘full’); % Branch active power flow variables
% Define the objective function
objective = sum(a2 .* Pg.^2 + a1 .* Pg + a0); % Minimize the total operating cost
% Define constraints
constraints = [];
% Node power balance constraints
for i = 1:nb
% Active power balance equation for node i
if ismember(i, gen_id)
record = find(gen_id == i);
constraints = [constraints, Pg(record) – Pd(i) – V(i) * sum(V * (Gbus(i, 🙂 * cos(theta(i) – theta) + Bbus(i, 🙂 * sin(theta(i) – theta))) == 0];
constraints = [constraints, Qg(record) – Qd(i) + V(i) * sum(V * (Gbus(i, 🙂 * sin(theta(i) – theta) – Bbus(i, 🙂 * cos(theta(i) – theta))) == 0];
else
constraints = [constraints, -Pd(i) – V(i) * sum(V * (Gbus(i, 🙂 * cos(theta(i) – theta) + Bbus(i, 🙂 * sin(theta(i) – theta))) == 0];
constraints = [constraints, -Qd(i) + V(i) * sum(V * (Gbus(i, 🙂 * sin(theta(i) – theta) – Bbus(i, 🙂 * cos(theta(i) – theta))) == 0];
end
end
% Generator output limit constraints
for i = 1:ng
% Minimum active power output constraint for generator i
constraints = [constraints, Pgmin(i) <= Pg(i)];
% Maximum active power output constraint for generator i
constraints = [constraints, Pg(i) <= Pgmax(i)];
% Minimum reactive power output constraint for generator i
constraints = [constraints, Qgmin(i) <= Qg(i)];
% Maximum reactive power output constraint for generator i
constraints = [constraints, Qg(i) <= Qgmax(i)];
end
% Bus voltage limit constraints
for i = 1:nb
% Minimum bus voltage magnitude constraint for bus i
constraints = [constraints, Vmin(i) <= V(i)];
% Maximum bus voltage magnitude constraint for bus i
constraints = [constraints, V(i) <= Vmax(i)];
end
% Branch power flow limit constraints
for i = 1:nl
% Calculate active power flow on branch i (assuming the "from" bus is the positive direction)
Pl(i) = V(fbus(i))^2 * Gbus(fbus(i), tbus(i)) – V(fbus(i)) * V(tbus(i)) * (Gbus(fbus(i), tbus(i)) * cos(theta(fbus(i) – theta(tbus(i))) + Bbus(fbus(i), tbus(i)) * sin(theta(fbus(i) – theta(tbus(i)))));
% Absolute value constraint for active power flow on branch i
constraints = [constraints, abs(Pl(i)) <= Plmax(i)];
end
% Solve the optimal power flow problem
result = solvesdp(constraints, objective); % Use YALMIP to solve the optimal power flow problem
% Output and analyze the results
if result.problem == 0 % If the solution is successful
disp(‘Optimal Power Flow Solved Successfully!’); % Display success message
disp(‘Total System Operating Cost:’); % Display total system operating cost
value(objective) % Display the objective function value
disp(‘Generator Active Power Output:’); % Display generator active power output
value(Pg) % Display the generator active power output variable values
disp(‘Generator Reactive Power Output:’); % Display generator reactive power output
value(Qg) % Display the generator reactive power output variable values
disp(‘Bus Voltage Magnitude:’); % Display bus voltage magnitude
value(V) % Display the bus voltage magnitude variable values
disp(‘Bus Voltage Angle (in degrees):’); % Display bus voltage angle (in degrees)
rad2deg(value(theta)) % Display bus voltage angle variable values and convert to degrees
disp(‘Branch Active Power Flow:’); % Display branch active power flow
value(Pl)
endI use MATLAB-YALMIP to write the optimal trend program based on the IEEEE14 system, but I can’t converge. Please help to see if the program is wrong?
code:
% Optimal Power Flow Example
% Based on the IEEE 14-bus standard case (i.e., case14 in Matpower)
% Implemented in MATLAB with YALMIP
clear; clc;
% Load power system data
mpc = loadcase(‘case14’); % You can also load other data files, such as case9, case30, etc.
% Define variables and parameters
nb = size(mpc.bus, 1); % Number of buses (nodes)
nl = size(mpc.branch, 1); % Number of branches (lines)
ng = size(mpc.gen, 1); % Number of generators
nr = nb – ng; % Number of load buses
% Extract node data
bus_id = mpc.bus(:, 1); % Bus IDs
bus_type = mpc.bus(:, 2); % Bus types (1-PQ, 2-PV, 3-Slack)
Pd = mpc.bus(:, 3); % Active power demand (MW)
Qd = mpc.bus(:, 4); % Reactive power demand (MVAr)
Vmin = mpc.bus(:, 13); % Minimum bus voltage (p.u.)
Vmax = mpc.bus(:, 12); % Maximum bus voltage (p.u.)
% Extract generator data
gen_id = mpc.gen(:, 1); % Bus IDs where generators are connected
Pgmin = mpc.gen(:, 10); % Minimum generator active power output (MW)
Pgmax = mpc.gen(:, 9); % Maximum generator active power output (MW)
Qgmin = mpc.gen(:, 5); % Minimum generator reactive power output (MVAr)
Qgmax = mpc.gen(:, 4); % Maximum generator reactive power output (MVAr)
a2 = mpc.gencost(:, 5); % Generator cost curve coefficient for quadratic term
a1 = mpc.gencost(:, 6); % Generator cost curve coefficient for linear term
a0 = mpc.gencost(:, 7); % Generator cost curve constant term
% Extract branch data
fbus = mpc.branch(:, 1); % "From" bus (starting node) of each branch
tbus = mpc.branch(:, 2); % "To" bus (ending node) of each branch
Plmax = ones(nl, 1) * 4000; % Maximum branch active power flow limit (MW)
% Plmax = mpc.branch(:, 6) .* mpc.branch(:, 3) * 1e-3; % Maximum branch active power flow limit (MW)
% Calculate the bus admittance matrix
[Ybus, ~, ~] = makeYbus(mpc); % Use a Matpower function to calculate the bus admittance matrix
Gbus = real(Ybus); % Real part of bus admittance matrix (conductance)
Bbus = imag(Ybus); % Imaginary part of bus admittance matrix (susceptance)
% Define YALMIP variables
Pg = sdpvar(ng, 1, ‘full’); % Generator active power output variables
Qg = sdpvar(ng, 1, ‘full’); % Generator reactive power output variables
V = sdpvar(nb, 1, ‘full’); % Bus voltage magnitude variables
theta = sdpvar(nb, 1, ‘full’); % Bus voltage angle variables
Pl = sdpvar(nl, 1, ‘full’); % Branch active power flow variables
% Define the objective function
objective = sum(a2 .* Pg.^2 + a1 .* Pg + a0); % Minimize the total operating cost
% Define constraints
constraints = [];
% Node power balance constraints
for i = 1:nb
% Active power balance equation for node i
if ismember(i, gen_id)
record = find(gen_id == i);
constraints = [constraints, Pg(record) – Pd(i) – V(i) * sum(V * (Gbus(i, 🙂 * cos(theta(i) – theta) + Bbus(i, 🙂 * sin(theta(i) – theta))) == 0];
constraints = [constraints, Qg(record) – Qd(i) + V(i) * sum(V * (Gbus(i, 🙂 * sin(theta(i) – theta) – Bbus(i, 🙂 * cos(theta(i) – theta))) == 0];
else
constraints = [constraints, -Pd(i) – V(i) * sum(V * (Gbus(i, 🙂 * cos(theta(i) – theta) + Bbus(i, 🙂 * sin(theta(i) – theta))) == 0];
constraints = [constraints, -Qd(i) + V(i) * sum(V * (Gbus(i, 🙂 * sin(theta(i) – theta) – Bbus(i, 🙂 * cos(theta(i) – theta))) == 0];
end
end
% Generator output limit constraints
for i = 1:ng
% Minimum active power output constraint for generator i
constraints = [constraints, Pgmin(i) <= Pg(i)];
% Maximum active power output constraint for generator i
constraints = [constraints, Pg(i) <= Pgmax(i)];
% Minimum reactive power output constraint for generator i
constraints = [constraints, Qgmin(i) <= Qg(i)];
% Maximum reactive power output constraint for generator i
constraints = [constraints, Qg(i) <= Qgmax(i)];
end
% Bus voltage limit constraints
for i = 1:nb
% Minimum bus voltage magnitude constraint for bus i
constraints = [constraints, Vmin(i) <= V(i)];
% Maximum bus voltage magnitude constraint for bus i
constraints = [constraints, V(i) <= Vmax(i)];
end
% Branch power flow limit constraints
for i = 1:nl
% Calculate active power flow on branch i (assuming the "from" bus is the positive direction)
Pl(i) = V(fbus(i))^2 * Gbus(fbus(i), tbus(i)) – V(fbus(i)) * V(tbus(i)) * (Gbus(fbus(i), tbus(i)) * cos(theta(fbus(i) – theta(tbus(i))) + Bbus(fbus(i), tbus(i)) * sin(theta(fbus(i) – theta(tbus(i)))));
% Absolute value constraint for active power flow on branch i
constraints = [constraints, abs(Pl(i)) <= Plmax(i)];
end
% Solve the optimal power flow problem
result = solvesdp(constraints, objective); % Use YALMIP to solve the optimal power flow problem
% Output and analyze the results
if result.problem == 0 % If the solution is successful
disp(‘Optimal Power Flow Solved Successfully!’); % Display success message
disp(‘Total System Operating Cost:’); % Display total system operating cost
value(objective) % Display the objective function value
disp(‘Generator Active Power Output:’); % Display generator active power output
value(Pg) % Display the generator active power output variable values
disp(‘Generator Reactive Power Output:’); % Display generator reactive power output
value(Qg) % Display the generator reactive power output variable values
disp(‘Bus Voltage Magnitude:’); % Display bus voltage magnitude
value(V) % Display the bus voltage magnitude variable values
disp(‘Bus Voltage Angle (in degrees):’); % Display bus voltage angle (in degrees)
rad2deg(value(theta)) % Display bus voltage angle variable values and convert to degrees
disp(‘Branch Active Power Flow:’); % Display branch active power flow
value(Pl)
end I use MATLAB-YALMIP to write the optimal trend program based on the IEEEE14 system, but I can’t converge. Please help to see if the program is wrong?
code:
% Optimal Power Flow Example
% Based on the IEEE 14-bus standard case (i.e., case14 in Matpower)
% Implemented in MATLAB with YALMIP
clear; clc;
% Load power system data
mpc = loadcase(‘case14’); % You can also load other data files, such as case9, case30, etc.
% Define variables and parameters
nb = size(mpc.bus, 1); % Number of buses (nodes)
nl = size(mpc.branch, 1); % Number of branches (lines)
ng = size(mpc.gen, 1); % Number of generators
nr = nb – ng; % Number of load buses
% Extract node data
bus_id = mpc.bus(:, 1); % Bus IDs
bus_type = mpc.bus(:, 2); % Bus types (1-PQ, 2-PV, 3-Slack)
Pd = mpc.bus(:, 3); % Active power demand (MW)
Qd = mpc.bus(:, 4); % Reactive power demand (MVAr)
Vmin = mpc.bus(:, 13); % Minimum bus voltage (p.u.)
Vmax = mpc.bus(:, 12); % Maximum bus voltage (p.u.)
% Extract generator data
gen_id = mpc.gen(:, 1); % Bus IDs where generators are connected
Pgmin = mpc.gen(:, 10); % Minimum generator active power output (MW)
Pgmax = mpc.gen(:, 9); % Maximum generator active power output (MW)
Qgmin = mpc.gen(:, 5); % Minimum generator reactive power output (MVAr)
Qgmax = mpc.gen(:, 4); % Maximum generator reactive power output (MVAr)
a2 = mpc.gencost(:, 5); % Generator cost curve coefficient for quadratic term
a1 = mpc.gencost(:, 6); % Generator cost curve coefficient for linear term
a0 = mpc.gencost(:, 7); % Generator cost curve constant term
% Extract branch data
fbus = mpc.branch(:, 1); % "From" bus (starting node) of each branch
tbus = mpc.branch(:, 2); % "To" bus (ending node) of each branch
Plmax = ones(nl, 1) * 4000; % Maximum branch active power flow limit (MW)
% Plmax = mpc.branch(:, 6) .* mpc.branch(:, 3) * 1e-3; % Maximum branch active power flow limit (MW)
% Calculate the bus admittance matrix
[Ybus, ~, ~] = makeYbus(mpc); % Use a Matpower function to calculate the bus admittance matrix
Gbus = real(Ybus); % Real part of bus admittance matrix (conductance)
Bbus = imag(Ybus); % Imaginary part of bus admittance matrix (susceptance)
% Define YALMIP variables
Pg = sdpvar(ng, 1, ‘full’); % Generator active power output variables
Qg = sdpvar(ng, 1, ‘full’); % Generator reactive power output variables
V = sdpvar(nb, 1, ‘full’); % Bus voltage magnitude variables
theta = sdpvar(nb, 1, ‘full’); % Bus voltage angle variables
Pl = sdpvar(nl, 1, ‘full’); % Branch active power flow variables
% Define the objective function
objective = sum(a2 .* Pg.^2 + a1 .* Pg + a0); % Minimize the total operating cost
% Define constraints
constraints = [];
% Node power balance constraints
for i = 1:nb
% Active power balance equation for node i
if ismember(i, gen_id)
record = find(gen_id == i);
constraints = [constraints, Pg(record) – Pd(i) – V(i) * sum(V * (Gbus(i, 🙂 * cos(theta(i) – theta) + Bbus(i, 🙂 * sin(theta(i) – theta))) == 0];
constraints = [constraints, Qg(record) – Qd(i) + V(i) * sum(V * (Gbus(i, 🙂 * sin(theta(i) – theta) – Bbus(i, 🙂 * cos(theta(i) – theta))) == 0];
else
constraints = [constraints, -Pd(i) – V(i) * sum(V * (Gbus(i, 🙂 * cos(theta(i) – theta) + Bbus(i, 🙂 * sin(theta(i) – theta))) == 0];
constraints = [constraints, -Qd(i) + V(i) * sum(V * (Gbus(i, 🙂 * sin(theta(i) – theta) – Bbus(i, 🙂 * cos(theta(i) – theta))) == 0];
end
end
% Generator output limit constraints
for i = 1:ng
% Minimum active power output constraint for generator i
constraints = [constraints, Pgmin(i) <= Pg(i)];
% Maximum active power output constraint for generator i
constraints = [constraints, Pg(i) <= Pgmax(i)];
% Minimum reactive power output constraint for generator i
constraints = [constraints, Qgmin(i) <= Qg(i)];
% Maximum reactive power output constraint for generator i
constraints = [constraints, Qg(i) <= Qgmax(i)];
end
% Bus voltage limit constraints
for i = 1:nb
% Minimum bus voltage magnitude constraint for bus i
constraints = [constraints, Vmin(i) <= V(i)];
% Maximum bus voltage magnitude constraint for bus i
constraints = [constraints, V(i) <= Vmax(i)];
end
% Branch power flow limit constraints
for i = 1:nl
% Calculate active power flow on branch i (assuming the "from" bus is the positive direction)
Pl(i) = V(fbus(i))^2 * Gbus(fbus(i), tbus(i)) – V(fbus(i)) * V(tbus(i)) * (Gbus(fbus(i), tbus(i)) * cos(theta(fbus(i) – theta(tbus(i))) + Bbus(fbus(i), tbus(i)) * sin(theta(fbus(i) – theta(tbus(i)))));
% Absolute value constraint for active power flow on branch i
constraints = [constraints, abs(Pl(i)) <= Plmax(i)];
end
% Solve the optimal power flow problem
result = solvesdp(constraints, objective); % Use YALMIP to solve the optimal power flow problem
% Output and analyze the results
if result.problem == 0 % If the solution is successful
disp(‘Optimal Power Flow Solved Successfully!’); % Display success message
disp(‘Total System Operating Cost:’); % Display total system operating cost
value(objective) % Display the objective function value
disp(‘Generator Active Power Output:’); % Display generator active power output
value(Pg) % Display the generator active power output variable values
disp(‘Generator Reactive Power Output:’); % Display generator reactive power output
value(Qg) % Display the generator reactive power output variable values
disp(‘Bus Voltage Magnitude:’); % Display bus voltage magnitude
value(V) % Display the bus voltage magnitude variable values
disp(‘Bus Voltage Angle (in degrees):’); % Display bus voltage angle (in degrees)
rad2deg(value(theta)) % Display bus voltage angle variable values and convert to degrees
disp(‘Branch Active Power Flow:’); % Display branch active power flow
value(Pl)
end yalmip, opf, matlab, optimal power flow MATLAB Answers — New Questions
