refactored code | #3, #6

KIT-IAI · Nov 19, 2020 · 8205ec7 · 8205ec7
1 parent 347db65
commit 8205ec7
Show file tree

Hide file tree

Showing 10 changed files with 148 additions and 123 deletions.
diff --git a/03_parser/auxfuns/modelfuns/generate_local_power_flow_problem.m b/03_parser/auxfuns/modelfuns/generate_local_power_flow_problem.m
@@ -21,49 +21,28 @@
 %%% always at the end of the bus numbering.
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     buses_core = mpc.(names.regions.global);
-    N_core = numel(buses_core);
+    Ncore = numel(buses_core);
     copy_buses_local = mpc.(names.copy_buses.local);
-    N_copy = numel(copy_buses_local);
-
-
-   %% optimal power flow cost + gradient + hessian
-   [cost, grad_cost, eq, eq_jac, ineq, ineq_jac, lagrangian_hessian, state] = build_local_cost_function(mpc, names, postfix);
+    Ncopy = numel(copy_buses_local);
 
+    %% preparation
+    [mpc_opf, om, local_buses_to_remove, mpopt] = prepare_case_file(mpc, names);
+    [constraint_function, lagrangian_hessian] = build_local_constraint_function(mpc_opf, om, mpopt);
+    %% cost function + cost gradient
+    [cost, grad_cost] = build_local_cost_function(om);
+    %% equalities + Jacobian
+    [eq, eq_jac] = build_local_equalities(constraint_function, local_buses_to_remove);
+    %% inequalities + Jacobian
+    [ineq, ineq_jac] = build_local_inequalities(constraint_function);
+    %% symbolic state
+    state = build_local_state(mpc_opf, names, postfix);
     %% initial conditions
     x0 = rand()
     %% lower and upper bounds
     % [lb, ub] = ...
 
-    %% dimensions
-    dims = [];
-    %% generate return values
-%     if strcmp(problem_type,'feasibility')
-%         cost = @(x) opf_p(x);
-%         grad_cost = @(x)gradient_costs(x);
-%         % grad_cost = @(x) zeros(4*N_core + 2*N_copy, 1);
-%         % TODO modify for OPF -> add second derivative of f
-%         Hessian = @(x, kappa, rho)jacobian_num(@(y)[Jac_pf(y); Jac_bus]'*kappa, x,  4*N_core + 2*N_copy, 4*N_core+ 2*N_copy);
-%         % cost = @(x) opf_p(x);
-%         ineq = @(x) [];
-%         eq = @(x)[ pf_p(x); pf_q(x); bus_specifications(x) ];
-%         pf = @(x)[ pf_p(x); pf_q(x) ];
-%         % TODO modify to gradient of cost (and later of h)
-%         Jac = Jac_g_ls;
-%         dims.eq = 4*N_core;
-%         dims.ineq = []; 
-%        % dims.ineq = length(mpc.gencost(:, 1));
-%     elseif strcmp(problem_type,'least-squares')
-%         g_ls    =  @(x)[pf_p(x); pf_q(x); bus_specifications(x)];
-%         grad_cost = @(x)(2*Jac_g_ls(x)'* g_ls(x));
-%         Hessian =  @(x,kappa, rho)(2*Jac_g_ls(x)'*Jac_g_ls(x));%@(x,kappa, rho)(2*Jac_g_ls(x)'*Jac_g_ls(x)); %@(x, kappa, rho)jacobian_num(@(y)(grad_cost(y)), x, 4*N_core + 2*N_copy, 4*N_core + 2*N_copy);
-%         cost = @(x)(g_ls(x)'*g_ls(x));
-%         ineq = @(x)[];
-%         eq = @(x)[];
-%         pf = @(x)[ pf_p(x); pf_q(x) ];
-%         Jac = @(x)[];
-%         dims.eq = [];
-%         dims.ineq = [];
-%     end
+    %% dimensions of state, equalities, inequalities
+    dims = build_local_dimensions(mpc_opf, ineq);
 end
 
 function entries = build_entries(N_core, N_copy, with_core)

diff --git a/06_opf_extension/build_local_constraint_function.m b/06_opf_extension/build_local_constraint_function.m
@@ -0,0 +1,7 @@
+function [constraint_function, Lxx] = build_local_constraint_function(mpc_opf, om, mpopt)
+    [Ybus, Yf, Yt] = makeYbus(mpc_opf);
+    il = find(mpc_opf.branch(:, 6) ~= 0 & mpc_opf.branch(:, 6) < 1e10);
+
+    constraint_function = @(x)opf_consfcn(x, om, Ybus, Yf(il,:), Yt(il,:), mpopt, il);
+    Lxx = @(x, lambda, cost_mult)opf_hessfcn(x, lambda, cost_mult, om, Ybus, Yf(il,:), Yt(il,:), mpopt, il);
+end
diff --git a/06_opf_extension/build_local_cost_function.m b/06_opf_extension/build_local_cost_function.m
@@ -1,73 +1,11 @@
-function [cost, grad, eq, eq_jac, ineq, ineq_jac, Lxx, state] = build_local_cost_function(mpc, names, postfix)
- [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
-            MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
-            QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
-
-
-    mpc = ext2int(mpc);
-    copy_buses_local = mpc.(names.copy_buses.local);
-    N_copy = numel(copy_buses_local);
-    copy_bus_types = get_bus_types(mpc.bus, copy_buses_local);
-
-    % turn off generators at copy nodes
-    for i = 1:length(copy_bus_types)
-        if copy_bus_types(i) == 2 || copy_bus_types(i) == 3
-            % get correct generator row in mpc.gen field 
-            gen_entry = find_generator_gen_entry(mpc, copy_buses_local(i));
-            % turn off corresponding generator
-            mpc.gen(gen_entry, GEN_STATUS) = 0;
-        end
-    end
-
-    % we changed the case file after it was switched to internal indexing
-    % we need to account for that
-    mpc.order.state = 'e';
-
-    % 2. generate cost function with opf_costfcn
-    [mpc_opf, mpopt] = opf_args(mpc); % only respect most simple opf formulation so far
-    mpc_opf = ext2int(mpc_opf);
-    om = opf_setup(mpc_opf, mpopt);
-
+function [cost, grad] = build_local_cost_function(om)
     f = @(x)opf_costfcn(x, om);
     cost = @(x)get_cost(x, f);
     grad = @(x)get_cost_gradient(x, f);
-%     hess = @(x)get_cost_hess(x, f);
-
-    [Ybus, Yf, Yt] = makeYbus(mpc_opf);
-    il = find(mpc_opf.branch(:, 6) ~= 0 & mpc_opf.branch(:, 6) < 1e10);
-
-    % equality constraints
-    gh_fcn = @(x)opf_consfcn(x, om, Ybus, Yf(il,:), Yt(il,:), mpopt, il);
-    eq = @(x)get_eq_cons(x, gh_fcn, copy_buses_local);
-    eq_jac = @(x)get_eq_cons_jacobian(x, gh_fcn);
-
-    % inequality constraints
-    ineq = @(x)get_ineq_cons(x, gh_fcn);
-    ineq_jac = @(x)get_ineq_cons_jacobian(x, gh_fcn);
-
-    % Hessian of Lagrangian
-    Lxx = @(x, lambda, cost_mult)opf_hessfcn(x, lambda, cost_mult, om, Ybus, Yf(il,:), Yt(il,:), mpopt, il);
-
-    buses_core = mpc.(names.regions.global);
-    N_core = numel(buses_core);
-
-    Ngen_on = sum(mpc.gen(:, GEN_STATUS));
-    [Vang_core, Vmag_core, Pg, Qg] = create_state_mp(postfix, N_core, Ngen_on);
-    [Vang_copy, Vmag_copy, ~, ~] = create_state_mp(strcat(postfix, '_copy'), N_copy, 0);
-
-    Vang = [Vang_core; Vang_copy];
-    Vmag = [Vmag_core; Vmag_copy];
-
-    state = stack_state(Vang, Vmag, Pg, Qg);
-end
-
-function types = get_bus_types(bus, buses)
-[PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
-        VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
-    types = bus(buses, BUS_TYPE);
+    % hess = @(x)get_cost_hess(x, f);
 end
 
-%% cost function
+% helper functions
 function f = get_cost(x, f_fcn)
     [f,~] = f_fcn(x);
 end
@@ -76,25 +14,6 @@
     [~, df] = f_fcn(x);
 end
 
-function d2f = get_cost_hess(x, f_fcn)
-    [~,~,d2f] = f_fcn(x);
-end
-%% equalities
-function g = get_eq_cons(x, gh_fcn, local_buses_to_remove)
-    [~,g,~,~] = gh_fcn(x);
-    % remove power flow equations for all copy buses
-    inds = [local_buses_to_remove; 2 * local_buses_to_remove]
-    g(inds) = [];
-end
-
-function dg = get_eq_cons_jacobian(x, gh_fcn)
-    [~,~,~,dg] = gh_fcn(x);
-end
-%% inequalities
-function h = get_ineq_cons(x, gh_fcn)
-    [h,~,~,~] = gh_fcn(x);
-end
-
-function dh = get_ineq_cons_jacobian(x, gh_fcn)
-    [~,~,dh,~] = gh_fcn(x);
-end
+% function d2f = get_cost_hess(x, f_fcn)
+%     [~,~,d2f] = f_fcn(x);
+% end
diff --git a/06_opf_extension/build_local_dimensions.m b/06_opf_extension/build_local_dimensions.m
@@ -0,0 +1,10 @@
+function dims = build_local_dimensions(mpc_opf, ineq)
+    Nbus = size(mpc_opf.bus, 1);
+    Ngen = size(mpc_opf.gen, 1);
+    dims.state = 2 * (Nbus + Ngen);
+    dims.eq = 2 * Nbus;
+    %% ToDo!
+    % verification missing!!!
+    x = rand(dims.state, 1);
+    dims.ineq = numel(ineq(x));
+end
diff --git a/06_opf_extension/build_local_equalities.m b/06_opf_extension/build_local_equalities.m
@@ -0,0 +1,16 @@
+function [eq, eq_jac] = build_local_equalities(constraint_function, local_buses_to_remove)
+    eq = @(x)get_eq_cons(x, constraint_function, copy_buses_local);
+    eq_jac = @(x)get_eq_cons_jacobian(x, constraint_function);
+end
+
+%% equalities
+function g = get_eq_cons(x, gh_fcn, local_buses_to_remove)
+    [~,g,~,~] = gh_fcn(x);
+    % remove power flow equations for all copy buses
+    inds = [local_buses_to_remove; 2 * local_buses_to_remove];
+    g(inds) = [];
+end
+
+function dg = get_eq_cons_jacobian(x, gh_fcn)
+    [~,~,~,dg] = gh_fcn(x);
+end
diff --git a/06_opf_extension/build_local_inequalities.m b/06_opf_extension/build_local_inequalities.m
@@ -0,0 +1,13 @@
+function [ineq, ineq_jac] = build_local_inequalities(constraint_function)
+    ineq = @(x)get_ineq_cons(x, constraint_function);
+    ineq_jac = @(x)get_ineq_cons_jacobian(x, constraint_function);
+end
+
+%% inequalities
+function h = get_ineq_cons(x, gh_fcn)
+    [h,~,~,~] = gh_fcn(x);
+end
+
+function dh = get_ineq_cons_jacobian(x, gh_fcn)
+    [~,~,dh,~] = gh_fcn(x);
+end
diff --git a/06_opf_extension/build_local_state.m b/06_opf_extension/build_local_state.m
@@ -0,0 +1,12 @@
+function state = build_local_state(mpc, names, postfix)
+    Ngen_on = size(mpc.gen, 1);
+    Ncopy = numel(mpc.(names.copy_buses.local));
+    Ncore = size(mpc.bus, 1) - Ncopy;
+    [Vang_core, Vmag_core, Pg, Qg] = create_state_mp(postfix, Ncore, Ngen_on);
+    [Vang_copy, Vmag_copy, ~, ~] = create_state_mp(strcat(postfix, '_copy'), Ncopy, 0);
+
+    Vang = [Vang_core; Vang_copy];
+    Vmag = [Vmag_core; Vmag_copy];
+
+    state = stack_state(Vang, Vmag, Pg, Qg);
+end
diff --git a/06_opf_extension/mpc_merge.m b/06_opf_extension/mpc_merge.m
@@ -21,6 +21,11 @@
 	8	2	300	98.61	0	0	1	1	0	230	1	1.1	0.9;
 	9	2	400	131.47	0	0	1	1	0	230	1	1.1	0.9;
 	10	1	0	0	0	0	1	1.1	0	230	1	1.1	0.9;
+	11	1	0	0	0	0	1	1	0	230	1	1.1	0.9;
+	12	1	300	98.61	0	0	1	1	0	230	1	1.1	0.9;
+	13	2	300	98.61	0	0	1	1	0	230	1	1.1	0.9;
+	14	2	400	131.47	0	0	1	1	0	230	1	1.1	0.9;
+	15	2	0	0	0	0	1	1	0	230	1	1.1	0.9;
 ];
 
 %% generator data
@@ -35,6 +40,9 @@
 	6	170	0	127.5	-127.5	1	100	1	170	0	0	0	0	0	0	0	0	0	0	0	0;
 	8	323.49	0	390	-390	1	100	1	520	0	0	0	0	0	0	0	0	0	0	0	0;
 	9	0	0	150	-150	1	100	1	200	0	0	0	0	0	0	0	0	0	0	0	0;
+	13	323.49	0	390	-390	1	100	1	520	0	0	0	0	0	0	0	0	0	0	0	0;
+	14	0	0	150	-150	1	100	1	200	0	0	0	0	0	0	0	0	0	0	0	0;
+	15	466.51	0	450	-450	1.1	100	1	600	0	0	0	0	0	0	0	0	0	0	0	0;
 ];
 
 %% branch data
@@ -53,6 +61,13 @@
 	8	9	0.00297	0.0297	0.00674	0	0	0	0	0	1	-360	360;
 	9	10	0.00297	0.0297	0.00674	240	240	240	0	0	1	-360	360;
 	1	10	0	0.00623	0	0	0	0	0.985	0	1	0	0;
+	11	12	0.00281	0.0281	0.00712	400	400	400	0	0	1	-360	360;
+	11	14	0.00304	0.0304	0.00658	0	0	0	0	0	1	-360	360;
+	11	15	0.00064	0.0064	0.03126	0	0	0	0	0	1	-360	360;
+	12	13	0.00108	0.0108	0.01852	0	0	0	0	0	1	-360	360;
+	13	14	0.00297	0.0297	0.00674	0	0	0	0	0	1	-360	360;
+	14	15	0.00297	0.0297	0.00674	240	240	240	0	0	1	-360	360;
+	9	11	0	0.00623	0	0	0	0	0.985	0	1	0	0;
 ];
 
 %%-----  OPF Data  -----%%
@@ -69,4 +84,7 @@
 	2	0	0	2	15	0;
 	2	0	0	2	30	0;
 	2	0	0	2	40	0;
+	2	0	0	2	30	0;
+	2	0	0	2	40	0;
+	2	0	0	2	10	0;
 ];
diff --git a/06_opf_extension/mpc_merge_split.m b/06_opf_extension/mpc_merge_split.m
@@ -21,6 +21,11 @@
 	8	2	300	98.61	0	0	1	1	0	230	1	1.1	0.9;
 	9	2	400	131.47	0	0	1	1	0	230	1	1.1	0.9;
 	10	1	0	0	0	0	1	1.1	0	230	1	1.1	0.9;
+	11	1	0	0	0	0	1	1	0	230	1	1.1	0.9;
+	12	1	300	98.61	0	0	1	1	0	230	1	1.1	0.9;
+	13	2	300	98.61	0	0	1	1	0	230	1	1.1	0.9;
+	14	2	400	131.47	0	0	1	1	0	230	1	1.1	0.9;
+	15	2	0	0	0	0	1	1	0	230	1	1.1	0.9;
 ];
 
 %% generator data
@@ -35,6 +40,9 @@
 	6	170	0	127.5	-127.5	1	100	1	170	0	0	0	0	0	0	0	0	0	0	0	0;
 	8	323.49	0	390	-390	1	100	1	520	0	0	0	0	0	0	0	0	0	0	0	0;
 	9	0	0	150	-150	1	100	1	200	0	0	0	0	0	0	0	0	0	0	0	0;
+	13	323.49	0	390	-390	1	100	1	520	0	0	0	0	0	0	0	0	0	0	0	0;
+	14	0	0	150	-150	1	100	1	200	0	0	0	0	0	0	0	0	0	0	0	0;
+	15	466.51	0	450	-450	1.1	100	1	600	0	0	0	0	0	0	0	0	0	0	0	0;
 ];
 
 %% branch data
@@ -53,6 +61,13 @@
 	8	9	0.00297	0.0297	0.00674	0	0	0	0	0	1	-360	360;
 	9	10	0.00297	0.0297	0.00674	240	240	240	0	0	1	-360	360;
 	1	10	0	0.00623	0	0	0	0	0.985	0	1	0	0;
+	11	12	0.00281	0.0281	0.00712	400	400	400	0	0	1	-360	360;
+	11	14	0.00304	0.0304	0.00658	0	0	0	0	0	1	-360	360;
+	11	15	0.00064	0.0064	0.03126	0	0	0	0	0	1	-360	360;
+	12	13	0.00108	0.0108	0.01852	0	0	0	0	0	1	-360	360;
+	13	14	0.00297	0.0297	0.00674	0	0	0	0	0	1	-360	360;
+	14	15	0.00297	0.0297	0.00674	240	240	240	0	0	1	-360	360;
+	9	11	0	0.00623	0	0	0	0	0.985	0	1	0	0;
 ];
 
 %%-----  OPF Data  -----%%
@@ -69,4 +84,7 @@
 	2	0	0	2	15	0;
 	2	0	0	2	30	0;
 	2	0	0	2	40	0;
+	2	0	0	2	30	0;
+	2	0	0	2	40	0;
+	2	0	0	2	10	0;
 ];
diff --git a/06_opf_extension/prepare_case_file.m b/06_opf_extension/prepare_case_file.m
@@ -0,0 +1,33 @@
+function [mpc_opf, om, copy_buses_local, mpopt] = prepare_case_file(mpc, names)
+    %%
+    [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
+            MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
+            QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
+
+    mpc = ext2int(mpc);
+    copy_buses_local = mpc.(names.copy_buses.local);
+    copy_bus_types = get_bus_types(mpc.bus, copy_buses_local);
+
+    %% turn off generators at copy nodes
+    for i = 1:length(copy_bus_types)
+        if copy_bus_types(i) == 2 || copy_bus_types(i) == 3
+            % get correct generator row in mpc.gen field 
+            gen_entry = find_generator_gen_entry(mpc, copy_buses_local(i));
+            % turn off corresponding generator
+            mpc.gen(gen_entry, GEN_STATUS) = 0;
+        end
+    end
+    %% we changed the case file after it was switched to internal indexing
+    % we need to account for that
+    mpc.order.state = 'e';
+    %% return values
+    [mpc_opf, mpopt] = opf_args(mpc); % only respect most simple opf formulation so far
+    mpc_opf = ext2int(mpc_opf);
+    om = opf_setup(mpc_opf, mpopt);
+end
+
+function types = get_bus_types(bus, buses)
+    [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
+        VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
+    types = bus(buses, BUS_TYPE);
+end