[wip] pair programming session with @bennerh, #3

KIT-IAI · Nov 18, 2020 · 18cd0fd · 18cd0fd
1 parent 613f06d
commit 18cd0fd
Show file tree

Hide file tree

Showing 5 changed files with 99 additions and 58 deletions.
diff --git a/03_parser/auxfuns/create_state_mp.m b/03_parser/auxfuns/create_state_mp.m
@@ -0,0 +1,23 @@
+function [vang, vmag, pg, qg] = create_state_mp(postfix, Nbus, Ngen)
+% create_state
+%
+%   `copy the declaration of the function in here (leave the ticks unchanged)`
+%
+%   _describe what the function does in the following line_
+%
+%   # Markdown formatting is supported
+%   Equations are possible to, e.g $a^2 + b^2 = c^2$.
+%   So are lists:
+%   - item 1
+%   - item 2
+%   ```matlab
+%   function y = square(x)
+%       x^2
+%   end
+%   ```
+%   See also: [run_case_file_splitter](run_case_file_splitter.md)
+    vang = sym(strcat('Va_', postfix, '_'), [Nbus 1], 'real');
+    vmag = sym(strcat('Vm_', postfix, '_'), [Nbus 1], 'real');
+    pg = sym(strcat('Pg_', postfix, '_'), [Ngen 1], 'real');
+    qg = sym(strcat('Qg_', postfix, '_'), [Ngen 1], 'real');
+end
diff --git a/03_parser/auxfuns/modelfuns/generate_local_power_flow_problem.m b/03_parser/auxfuns/modelfuns/generate_local_power_flow_problem.m
@@ -1,4 +1,4 @@
-function [cost, ineq, eq, x0, grad_cost, eq_jac, ineq_jac, Hessian, state, dims] = generate_local_power_flow_problem(mpc, names, postfix, problem_type)
+function [cost, ineq, eq, x0, grad_cost, eq_jac, ineq_jac, lagrangian_hessian, state, dims] = generate_local_power_flow_problem(mpc, names, postfix, problem_type)
 % generate_local_power_flow_problem
 %
 %   `copy the declaration of the function in here (leave the ticks unchanged)`
@@ -25,50 +25,45 @@
     copy_buses_local = mpc.(names.copy_buses.local);
     N_copy = numel(copy_buses_local);
 
-    [Vang_core, Vmag_core, Pnet_core, Qnet_core] = create_state(postfix, N_core);
-    [Vang_copy, Vmag_copy, Pnet_copy, Qnet_copy] = create_state(strcat(postfix, '_copy'), N_copy);
 
-    Vang = [Vang_core; Vang_copy];
-    Vmag = [Vmag_core; Vmag_copy];
-    Pnet = Pnet_core;
-    Qnet = Qnet_core;
-
-    state = stack_state(Vang, Vmag, Pnet, Qnet);
    %% optimal power flow cost + gradient + hessian
-   [cost, grad_cost, hess_cost, eq, eq_jac, ineq, ineq_jac] = build_local_cost_function(mpc, names);
+   [cost, grad_cost, eq, eq_jac, ineq, ineq_jac, lagrangian_hessian, state] = build_local_cost_function(mpc, names, postfix);
 
     %% initial conditions
-    % x0 = ...
+    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
+%     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
 end
 
 function entries = build_entries(N_core, N_copy, with_core)

diff --git a/03_parser/generate_distributed_problem.m b/03_parser/generate_distributed_problem.m
@@ -18,16 +18,18 @@
 %   See also: [run_case_file_splitter](run_case_file_splitter.md)
     % extract Data from casefile
     [N_regions, N_buses_in_regions, N_copy_buses_in_regions, ~] = get_relevant_information(mpc, names);
-    [costs, inequalities, equalities, states, xx0, pfs, bus_specs, Jacs, grads, Hessians, dims] = deal(cell(N_regions,1));
+    [costs,  inequalities, equalities, xx0, grads, Jacs, Hessians, states, dims] = deal(cell(N_regions,1));
     connection_table = mpc.(names.consensus);
     % set up the Ai's
     consensus_matrices = create_consensus_matrices(connection_table, N_buses_in_regions, N_copy_buses_in_regions);
     % create local power flow problems
     fprintf('\n\n');
     for i = 1:N_regions
         fprintf('Creating power flow problem for system %i...', i);
-        [cost, inequality, equality, x0, pf, bus_spec, Jac, grad, Hessian, state, dim] = generate_local_power_flow_problem(mpc.(names.split){i}, names, num2str(i), problem_type);
-        [costs{i}, inequalities{i}, equalities{i}, xx0{i}, pfs{i}, bus_specs{i}, states{i}, Jacs{i}, grads{i}, Hessians{i}, dims{i}] = deal(cost, inequality, equality, x0, pf, bus_spec, state, Jac, grad, Hessian, dim);
+        [cost, inequality, equality, x0, grad, eq_jac, ineq_jac, Hessian, state, dim] = generate_local_power_flow_problem(mpc.(names.split){i}, names, num2str(i), problem_type);
+        % combine Jacobians of inequalities and equalities in single Jacobian
+        Jac = @(x)[eq_jac(x); ineq_jac(x)];
+        [costs{i},  inequalities{i}, equalities{i}, xx0{i}, grads{i}, Jacs{i}, Hessians{i}, states{i}, dims{i}] = deal(cost, inequality, equality, x0, grad, Jac, Hessian, state, dim);
         fprintf('done.\n')
     end
     %% generate output
@@ -44,8 +46,6 @@
     problem.zz0 = xx0;
     problem.AA  = consensus_matrices;
 
-    problem.pf = pfs;
-    problem.bus_specs = bus_specs;
     problem.state = states;
 end
 

diff --git a/06_opf_extension/build_local_cost_function.m b/06_opf_extension/build_local_cost_function.m
@@ -1,11 +1,15 @@
-function [cost, grad, hess, eq, eq_jac, ineq, ineq_jac] = build_local_cost_function(mpc, names)
+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;
 
-    copy_buses_local = mpc.(names.copy_buses.local);    
+
+    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 
@@ -15,32 +19,46 @@
         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, mpopt);
+    mpc_opf = ext2int(mpc_opf);
     om = opf_setup(mpc_opf, mpopt);
 
     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);
-
+%     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);
 
-    %% WARNING!!!!!!!!
-    % We need to delete the power flow equations for the copy buses!!!!!!!
-    %%
-
     % equality constraints
     gh_fcn = @(x)opf_consfcn(x, om, Ybus, Yf(il,:), Yt(il,:), mpopt, il);
-    eq = @(x)get_eq_cons(x, gh_fcn);
+    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)
@@ -62,8 +80,11 @@
     [~,~,d2f] = f_fcn(x);
 end
 %% equalities
-function g = get_eq_cons(x, gh_fcn)
+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)

diff --git a/06_opf_extension/opf_testfile.m b/06_opf_extension/opf_testfile.m
@@ -14,10 +14,12 @@
 mpc.fields_to_merge = {'bus', 'gen', 'branch', 'gencost'};
 %mpc.fields_to_merge = {'bus', 'gen', 'branch'};
 mpc.trans = loadcase('case5');
-mpc.dist = { loadcase('case5')  
+mpc.dist = { loadcase('case5');
+             loadcase('case5')
                     };
 
-mpc.connection_array = [ 1 2 1 5];
+mpc.connection_array = [ 1 2 1 5;
+                         2 3 4 1];
 
 % compatibility tests
 % Opf data is provided ??
@@ -79,8 +81,8 @@
 mpc_split = run_case_file_splitter(mpc_merge, conn, names);
 
 %% run opf of merged file
-runopf(mpc_merge);
-runpf(mpc_merge);
+% runopf(mpc_merge);
+% runpf(mpc_merge);
 %% setup distributed opf
 % start values from pf