Refactor dc sensivity analysis

src/main/java/com/powsybl/openloadflow/dc/DcLoadFlowEngine.java

@@ -10,10 +10,7 @@
 import com.powsybl.math.matrix.MatrixFactory;
 import com.powsybl.openloadflow.dc.equations.DcEquationSystem;
 import com.powsybl.openloadflow.dc.equations.DcEquationSystemCreationParameters;
-import com.powsybl.openloadflow.equations.EquationSystem;
-import com.powsybl.openloadflow.equations.JacobianMatrix;
-import com.powsybl.openloadflow.equations.UniformValueVoltageInitializer;
-import com.powsybl.openloadflow.equations.VariableSet;
+import com.powsybl.openloadflow.equations.*;
 import com.powsybl.openloadflow.network.FirstSlackBusSelector;
 import com.powsybl.openloadflow.network.LfBus;
 import com.powsybl.openloadflow.network.LfNetwork;
@@ -22,10 +19,7 @@
 import org.slf4j.Logger;
 import org.slf4j.LoggerFactory;
 
-import java.util.Arrays;
-import java.util.Collections;
-import java.util.List;
-import java.util.Objects;
+import java.util.*;
 
 /**
  * @author Geoffroy Jamgotchian <geoffroy.jamgotchian at rte-france.com>
@@ -38,6 +32,8 @@ public class DcLoadFlowEngine {
 
     private final DcLoadFlowParameters parameters;
 
+    private double[] targetVector;
+
     public DcLoadFlowEngine(LfNetwork network, MatrixFactory matrixFactory) {
         this.networks = Collections.singletonList(network);
         parameters = new DcLoadFlowParameters(new FirstSlackBusSelector(), matrixFactory);
@@ -53,52 +49,88 @@ public DcLoadFlowEngine(List<LfNetwork> networks, DcLoadFlowParameters parameter
         this.parameters = Objects.requireNonNull(parameters);
     }
 
-    private void distributeSlack(LfNetwork network) {
-        double mismatch = network.getActivePowerMismatch();
+    private void distributeSlack(Collection<LfBus> buses) {
+        double mismatch = getActivePowerMismatch(buses);
         ActivePowerDistribution activePowerDistribution = ActivePowerDistribution.create(parameters.getBalanceType(), false);
-        activePowerDistribution.run(network, mismatch);
+        activePowerDistribution.run(buses, mismatch);
+    }
+
+    private static double getActivePowerMismatch(Collection<LfBus> buses) {
+        double mismatch = 0;
+        for (LfBus b : buses) {
+            mismatch += b.getGenerationTargetP() - b.getLoadTargetP();
+        }
+        return -mismatch;
     }
 
     public DcLoadFlowResult run() {
         // only process main (largest) connected component
         LfNetwork network = networks.get(0);
 
-        if (parameters.isDistributedSlack()) {
-            distributeSlack(network);
-        }
-
         DcEquationSystemCreationParameters creationParameters = new DcEquationSystemCreationParameters(parameters.isUpdateFlows(), false, parameters.isForcePhaseControlOffAndAddAngle1Var(), parameters.isUseTransformerRatio());
         EquationSystem equationSystem = DcEquationSystem.create(network, new VariableSet(), creationParameters);
 
+        LoadFlowResult.ComponentResult.Status status = LoadFlowResult.ComponentResult.Status.FAILED;
+        try (JacobianMatrix j = new JacobianMatrix(equationSystem, parameters.getMatrixFactory())) {
+
+            status = run(equationSystem, j, Collections.emptyList());
+        } catch (Exception e) {
+            LOGGER.error("Failed to solve linear system for DC load flow", e);
+        }
+
+        return new DcLoadFlowResult(network, getActivePowerMismatch(network.getBuses()), status);
+    }
+
+    public LoadFlowResult.ComponentResult.Status run(EquationSystem equationSystem, JacobianMatrix j, Collection<LfBus> disabledBuses) {
+
         double[] x = equationSystem.createStateVector(new UniformValueVoltageInitializer());
 
-        equationSystem.updateEquations(x);
+        // only process main (largest) connected component
+        LfNetwork network = networks.get(0);
 
-        double[] targets = equationSystem.createTargetVector();
+        Collection<LfBus> remainingBuses = new HashSet<>(network.getBuses());
+        remainingBuses.removeAll(disabledBuses);
 
-        try (JacobianMatrix j = new JacobianMatrix(equationSystem, parameters.getMatrixFactory())) {
-            double[] dx = Arrays.copyOf(targets, targets.length);
+        if (parameters.isDistributedSlack()) {
+            distributeSlack(remainingBuses);
+        }
 
-            LoadFlowResult.ComponentResult.Status status;
-            try {
-                j.solveTransposed(dx);
-                status = LoadFlowResult.ComponentResult.Status.CONVERGED;
-            } catch (Exception e) {
-                status = LoadFlowResult.ComponentResult.Status.FAILED;
-                LOGGER.error("Failed to solve linear system for DC load flow", e);
-            }
+        equationSystem.updateEquations(x);
 
-            equationSystem.updateEquations(dx);
-            equationSystem.updateNetwork(dx);
+        this.targetVector = equationSystem.createTargetVector();
 
-            // set all calculated voltages to NaN
-            for (LfBus bus : network.getBuses()) {
-                bus.setV(Double.NaN);
-            }
+        if (!disabledBuses.isEmpty()) {
+            // set buses injections and transformers to 0
+            disabledBuses.stream()
+                .map(lfBus -> equationSystem.getEquation(lfBus.getNum(), EquationType.BUS_P))
+                .filter(Optional::isPresent)
+                .map(Optional::get)
+                .map(Equation::getColumn)
+                .forEach(column -> targetVector[column] = 0);
+        }
+
+        LoadFlowResult.ComponentResult.Status status;
+        try {
+            j.solveTransposed(targetVector);
+            status = LoadFlowResult.ComponentResult.Status.CONVERGED;
+        } catch (Exception e) {
+            status = LoadFlowResult.ComponentResult.Status.FAILED;
+            LOGGER.error("Failed to solve linear system for DC load flow", e);
+        }
 
-            LOGGER.info("Dc loadflow complete (status={})", status);
+        equationSystem.updateEquations(targetVector);
+        equationSystem.updateNetwork(targetVector);
 
-            return new DcLoadFlowResult(network, network.getActivePowerMismatch(), status);
+        // set all calculated voltages to NaN
+        for (LfBus bus : network.getBuses()) {
+            bus.setV(Double.NaN);
         }
+
+        LOGGER.info("Dc loadflow complete (status={})", status);
+        return status;
+    }
+
+    public double[] getTargetVector() {
+        return targetVector;
     }
 }

src/main/java/com/powsybl/openloadflow/graph/EvenShiloachGraphDecrementalConnectivity.java

@@ -16,6 +16,7 @@
 import org.slf4j.LoggerFactory;
 
 import java.util.*;
+import java.util.stream.Collectors;
 
 /**
  * Implementing the Even-Shiloach algorithm (see https://dl.acm.org/doi/10.1145/322234.322235)
@@ -131,6 +132,7 @@ public void reset() {
 
     @Override
     public int getComponentNumber(V vertex) {
+        checkVertex(vertex);
         lazyComputeConnectivity();
         updateVertexMapCache();
         return vertexToConnectedComponent.get(vertex);
@@ -142,6 +144,31 @@ public List<Set<V>> getSmallComponents() {
         return newConnectedComponents;
     }
 
+    @Override
+    public Set<V> getConnectedComponent(V vertex) {
+        checkVertex(vertex);
+        lazyComputeConnectivity();
+        updateVertexMapCache();
+        int cn = vertexToConnectedComponent.get(vertex);
+        return cn == 0 ? getMainConnectedComponent() : newConnectedComponents.get(cn - 1);
+    }
+
+    private Set<V> getMainConnectedComponent() {
+        return vertices.stream().filter(v -> newConnectedComponents.stream().noneMatch(cc -> cc.contains(v))).collect(Collectors.toSet());
+    }
+
+    @Override
+    public Set<V> getNonConnectedVertices(V vertex) {
+        checkVertex(vertex);
+        lazyComputeConnectivity();
+        List<Set<V>> nonConnectedComponents = new ArrayList<>(newConnectedComponents);
+        newConnectedComponents.stream().filter(c -> c.contains(vertex)).findFirst().ifPresent(c -> {
+            nonConnectedComponents.remove(c);
+            nonConnectedComponents.add(getMainConnectedComponent());
+        });
+        return nonConnectedComponents.stream().flatMap(Collection::stream).collect(Collectors.toSet());
+    }
+
     private void lazyComputeConnectivity() {
         if (init && unprocessedCutEdges.isEmpty()) {
             return;
@@ -218,6 +245,12 @@ private void updateVertexMapCache() {
         }
     }
 
+    private void checkVertex(V vertex) {
+        if (!graph.containsVertex(vertex)) {
+            throw new AssertionError("given vertex " + vertex + " is not in the graph");
+        }
+    }
+
     private interface GraphProcess {
         void next();
 

src/main/java/com/powsybl/openloadflow/graph/GraphDecrementalConnectivity.java

@@ -44,4 +44,8 @@ public interface GraphDecrementalConnectivity<V> {
      * @return the collection of small connected components
      */
     Collection<Set<V>> getSmallComponents();
+
+    Set<V> getConnectedComponent(V vertex);
+
+    Set<V> getNonConnectedVertices(V vertex);
 }

src/main/java/com/powsybl/openloadflow/graph/MinimumSpanningTreeGraphDecrementalConnectivity.java

@@ -91,33 +91,56 @@ private void invalidateMst() {
 
     @Override
     public int getComponentNumber(V vertex) {
-        if (this.mst == null) {
-            this.mst = new KruskalMinimumSpanningTrees().getSpanningTree();
-        }
-        if (this.parentMap == null) {
-            this.parentMap = mst.forest.getParentMap();
-            this.sortedRoots = mst.forest.getSortedRoots();
-        }
+        checkVertex(vertex);
+        lazyCompute();
         return sortedRoots.indexOf(parentMap.get(vertex));
     }
 
     @Override
     public List<Set<V>> getSmallComponents() {
+        List<Set<V>> components = getConnectedComponents();
+        return components.subList(1, components.size());
+    }
+
+    private List<Set<V>> getConnectedComponents() {
+        lazyCompute();
+        List<Set<V>> components = new ArrayList<>();
+        for (V root : sortedRoots) {
+            Set<V> set = parentMap.entrySet().stream()
+                .filter(e -> e.getValue() == root)
+                .map(Map.Entry::getKey).collect(Collectors.toSet());
+            components.add(set);
+        }
+        return components;
+    }
+
+    @Override
+    public Set<V> getConnectedComponent(V vertex) {
+        checkVertex(vertex);
+        return getConnectedComponents().get(getComponentNumber(vertex));
+    }
+
+    @Override
+    public Set<V> getNonConnectedVertices(V vertex) {
+        checkVertex(vertex);
+        return getConnectedComponents().stream().filter(component -> !component.contains(vertex))
+            .flatMap(Collection::stream).collect(Collectors.toSet());
+    }
+
+    private void lazyCompute() {
         if (this.mst == null) {
             this.mst = new KruskalMinimumSpanningTrees().getSpanningTree();
         }
         if (this.parentMap == null) {
             this.parentMap = mst.forest.getParentMap();
             this.sortedRoots = mst.forest.getSortedRoots();
         }
-        List<Set<V>> components = new ArrayList<>();
-        for (V root : sortedRoots) {
-            Set<V> set = parentMap.entrySet().stream()
-                .filter(e -> e.getValue() == root)
-                .map(Map.Entry::getKey).collect(Collectors.toSet());
-            components.add(set);
+    }
+
+    private void checkVertex(V vertex) {
+        if (!graph.containsVertex(vertex)) {
+            throw new AssertionError("given vertex " + vertex + " is not in the graph");
         }
-        return components.subList(1, components.size());
     }
 
     class KruskalMinimumSpanningTrees implements SpanningTreeAlgorithm<Object> {

src/main/java/com/powsybl/openloadflow/graph/NaiveGraphDecrementalConnectivity.java

@@ -92,6 +92,7 @@ public void reset() {
 
     @Override
     public int getComponentNumber(V vertex) {
+        checkVertex(vertex);
         updateComponents();
         return components[numGetter.applyAsInt(vertex)];
     }
@@ -101,4 +102,25 @@ public Collection<Set<V>> getSmallComponents() {
         updateComponents();
         return componentSets.subList(1, componentSets.size());
     }
+
+    @Override
+    public Set<V> getConnectedComponent(V vertex) {
+        checkVertex(vertex);
+        updateComponents();
+        return componentSets.get(components[numGetter.applyAsInt(vertex)]);
+    }
+
+    @Override
+    public Set<V> getNonConnectedVertices(V vertex) {
+        checkVertex(vertex);
+        updateComponents();
+        return componentSets.stream().filter(component -> !component.contains(vertex))
+            .flatMap(Collection::stream).collect(Collectors.toSet());
+    }
+
+    private void checkVertex(V vertex) {
+        if (!graph.containsVertex(vertex)) {
+            throw new AssertionError("given vertex " + vertex + " is not in the graph");
+        }
+    }
 }

src/main/java/com/powsybl/openloadflow/network/LfNetwork.java

@@ -369,25 +369,6 @@ public void logBalance() {
                 num, activeGeneration, activeLoad, reactiveGeneration, reactiveLoad);
     }
 
-    public double getActivePowerMismatch() {
-        double mismatch = 0;
-        for (LfBus b : busesById.values()) {
-            mismatch += b.getGenerationTargetP() - b.getLoadTargetP();
-        }
-        return -mismatch;
-    }
-
-    public double getActivePowerMismatchInMainComponent(GraphDecrementalConnectivity<LfBus> connectivity) {
-        double mismatch = 0;
-        int mainComponent = connectivity.getComponentNumber(getSlackBus());
-        for (LfBus b : busesById.values()) {
-            if (connectivity.getComponentNumber(b) == mainComponent) {
-                mismatch += b.getGenerationTargetP() - b.getLoadTargetP();
-            }
-        }
-        return -mismatch;
-    }
-
     private static void fix(LfNetwork network, boolean minImpedance) {
         if (minImpedance) {
             for (LfBranch branch : network.getBranches()) {

src/main/java/com/powsybl/openloadflow/network/util/ActivePowerDistribution.java

@@ -7,8 +7,10 @@
 package com.powsybl.openloadflow.network.util;
 
 import com.powsybl.loadflow.LoadFlowParameters;
+import com.powsybl.openloadflow.network.LfBus;
 import com.powsybl.openloadflow.network.LfNetwork;
 
+import java.util.Collection;
 import java.util.List;
 import java.util.Objects;
 
@@ -26,7 +28,7 @@ public interface Step {
 
         String getElementType();
 
-        List<ParticipatingElement> getParticipatingElements(LfNetwork network);
+        List<ParticipatingElement> getParticipatingElements(Collection<LfBus> buses);
 
         double run(List<ParticipatingElement> participatingElements, int iteration, double remainingMismatch);
     }
@@ -62,7 +64,11 @@ public String getElementType() {
     }
 
     public Result run(LfNetwork network, double activePowerMismatch) {
-        List<ParticipatingElement> participatingElements = step.getParticipatingElements(network);
+        return run(network.getBuses(), activePowerMismatch);
+    }
+
+    public Result run(Collection<LfBus> buses, double activePowerMismatch) {
+        List<ParticipatingElement> participatingElements = step.getParticipatingElements(buses);
 
         int iteration = 0;
         double remainingMismatch = activePowerMismatch;