pair programming w/ @bennerh

- closes #7
- pertains #3
- identified and openend #11
- breaks with #10

00_use-case/getting_started_opf.m

@@ -90,7 +90,7 @@
 % generate distributed problem
 problem = generate_distributed_opf_for_aladin(mpc_split, names, 'feasibility');
 problem.solver = 'fmincon';
-[xval, xval_stacked] = validate_distributed_problem_formulation(problem, mpc_split, names);
+% [xval, xval_stacked] = validate_distributed_problem_formulation(problem, mpc_split, names);
 
 % solve distributed ALADIN
 

03_parser/auxfuns/create_consensus_matrices_opf.m

@@ -0,0 +1,81 @@
+function A = create_consensus_matrices_opf(tab, number_of_buses_in_region, number_of_generators_in_region)
+% create_consensus_matrices_opf
+%
+%   `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)
+    assert(istable(tab), 'expecting tab to be a table.');
+    assert(mod(height(tab), 2) == 0, 'inconsistent number of consensus restrictions.')
+    assert(numel(number_of_buses_in_region) == numel(number_of_generators_in_region), 'inconsistent dimensions');
+
+    % consensus for voltage angles
+    Aang = build_consensus_matrix_core(tab, number_of_buses_in_region, number_of_generators_in_region, 'ang');
+    % consensus for voltage magnitudes
+    Amag = build_consensus_matrix_core(tab, number_of_buses_in_region, number_of_generators_in_region, 'mag');
+    % stack together and check for correct dimensions
+    A = stack_and_check(Aang, Amag, height(tab), number_of_buses_in_region, number_of_generators_in_region);
+end
+
+function A = stack_and_check(Aang, Amag, Nconsensus, number_of_buses_in_region, number_of_copy_buses_in_region)
+    assert(numel(Aang) == numel(Amag), 'inconsistent dimensions.');
+    A = Aang;
+    Nrows = 2*Nconsensus;
+    for i = 1:numel(Amag)
+        A{i} = [ Aang{i}; Amag{i} ];
+%         Ncopy = number_of_copy_buses_in_region(i);
+%         Ncore = number_of_buses_in_region(i) - Ncopy;
+%         Ncols = 4*Ncore + 2*Ncopy;
+%         
+%         assert(prod(size(A{i}) == [Nrows, Ncols]) == 1, 'inconsistent dimensions for consensus matrix in region %i', i);
+    end
+end
+
+function A = build_consensus_matrix_core(tab, number_of_buses_in_region, number_of_generators_in_region, kind)
+    kind = lower(kind);
+    Nconsensus = height(tab);
+    Nregions = numel(number_of_buses_in_region);
+    A = cell(Nregions, 1);
+    for i = 1:Nregions
+        A{i} = sparse(Nconsensus, 2*number_of_buses_in_region(i) + 2*number_of_generators_in_region(i));
+    end
+
+    for i = 1:Nconsensus
+        orig_sys = tab.orig_sys(i);
+        copy_sys = tab.copy_sys(i);
+        orig_bus = tab.orig_bus_local(i);
+        copy_bus = tab.copy_bus_local(i);
+
+        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+        %%% the following code (implictly) assumes that the state in each
+        %%% region is stacked according to x = [vang, vmag, pnet, qnet],
+        %%% hence some shifting might be necessary, depending on whether we
+        %%% build consensus for the voltage angle or voltage magnitude.
+        if kind == 'ang'
+            Norig = 0;
+            Ncopy = 0;
+        elseif kind == 'mag'
+            Norig = number_of_buses_in_region(orig_sys);
+            Ncopy = number_of_buses_in_region(copy_sys);
+        else
+            error('kind %s is not supported', kind);
+        end
+        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+        orig_bus_entry = orig_bus + Norig;
+        copy_bus_entry = copy_bus + Ncopy;
+
+        A{orig_sys}(i, orig_bus_entry) = 1;
+        A{copy_sys}(i, copy_bus_entry) = -1;
+    end
+end

03_parser/auxfuns/modelfuns/build_local_opf.m

@@ -36,5 +36,5 @@
     %% lower and upper bounds
     [lb, ub] = build_local_bounds(om);
     %% dimensions of state, equalities, inequalities
-    dims = build_local_dimensions(mpc_opf, ineq);
+    dims = build_local_dimensions(mpc_opf, eq, ineq);
 end

03_parser/generate_distributed_opf.m

@@ -21,17 +21,20 @@
     [costs,  inequalities, equalities, xx0, grads, Jacs, Hessians, states, dims, lbs, ubs] = 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, grad, eq_jac, ineq_jac, Hessian, state, dim, lb, ub] = build_local_opf(mpc.(names.split){i}, names, num2str(i));
         % combine Jacobians of inequalities and equalities in single Jacobian
-        Jac = @(x)[eq_jac(x); ineq_jac(x)];
+        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}, lbs{i}, ubs{i}] = deal(cost, inequality, equality, x0, grad, Jac, Hessian, state, dim, lb, ub);
         fprintf('done.\n')
     end
+
+    N_generators_in_regions = struct_for_N_generators(dims);
+    consensus_matrices = create_consensus_matrices_opf(connection_table, N_buses_in_regions, N_generators_in_regions);
     %% generate output for Aladin
     problem.locFuns.ffi = costs;
     problem.locFuns.ggi = equalities;
@@ -51,6 +54,14 @@
     problem.llbx = lbs;
     problem.uubx = ubs;
 end
+
+function N_generators = struct_for_N_generators(dims)
+    n = numel(dims);
+    N_generators = zeros(n, 1);
+    for i = 1:n
+        N_generators(i) = dims{i}.n.gen;
+    end
+end
 
 
 

06_opf_extension/build_local_constraint_function.m

@@ -3,5 +3,8 @@
     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);
+    % cost multiplier is set to one
+    cost_mult = 1;
+    % aladin needs the rho value for Lxx
+    Lxx = @(x, lambda, rho)opf_hessfcn(x, lambda, cost_mult, om, Ybus, Yf(il,:), Yt(il,:), mpopt, il);
 end

06_opf_extension/build_local_dimensions.m

@@ -1,10 +1,15 @@
-function dims = build_local_dimensions(mpc_opf, ineq)
+function dims = build_local_dimensions(mpc_opf, eq, ineq)
+    % this file assumes that all generators at the copy busses have been turned off
+    % i.e. prepare_case_file() was called already!
     Nbus = size(mpc_opf.bus, 1);
     Ngen = size(mpc_opf.gen, 1);
     dims.state = 2 * (Nbus + Ngen);
     dims.eq = 2 * Nbus;
+    dims.n.bus = Nbus;
+    dims.n.gen = Ngen;
     %% ToDo!
     % verification missing!!!
     x = rand(dims.state, 1);
     dims.ineq = numel(ineq(x));
+    dims.eq = numel(eq(x));
 end

06_opf_extension/build_local_equalities.m

@@ -1,16 +1,17 @@
 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);
+    inds = [local_buses_to_remove; 2 * local_buses_to_remove];
+    eq = @(x)get_eq_cons(x, constraint_function, inds);
+    eq_jac = @(x)get_eq_cons_jacobian(x, constraint_function, inds);
 end
 
 %% equalities
-function g = get_eq_cons(x, gh_fcn, local_buses_to_remove)
+function g = get_eq_cons(x, gh_fcn, inds)
     [~,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)
+function dg = get_eq_cons_jacobian(x, gh_fcn, inds)
     [~,~,~,dg] = gh_fcn(x);
+    dg(:, inds) = [];
 end