Skip to content

Commit

Permalink
Refactor equation vector (#575)
Browse files Browse the repository at this point in the history
Signed-off-by: Geoffroy Jamgotchian <[email protected]>
  • Loading branch information
geofjamg authored Jul 2, 2022
1 parent 9153bfb commit e6a160a
Show file tree
Hide file tree
Showing 9 changed files with 255 additions and 142 deletions.
50 changes: 33 additions & 17 deletions src/main/java/com/powsybl/openloadflow/ac/nr/NewtonRaphson.java
Original file line number Diff line number Diff line change
Expand Up @@ -13,10 +13,15 @@
import com.powsybl.openloadflow.network.LfElement;
import com.powsybl.openloadflow.network.LfNetwork;
import com.powsybl.openloadflow.network.util.VoltageInitializer;
import org.apache.commons.lang3.tuple.Pair;
import org.slf4j.Logger;
import org.slf4j.LoggerFactory;

import java.util.Comparator;
import java.util.List;
import java.util.Map;
import java.util.Objects;
import java.util.stream.Collectors;

/**
* @author Geoffroy Jamgotchian <geoffroy.jamgotchian at rte-france.com>
Expand All @@ -37,38 +42,52 @@ public class NewtonRaphson {

private final TargetVector<AcVariableType, AcEquationType> targetVector;

private final EquationVector<AcVariableType, AcEquationType> equationVector;

public NewtonRaphson(LfNetwork network, NewtonRaphsonParameters parameters,
EquationSystem<AcVariableType, AcEquationType> equationSystem, JacobianMatrix<AcVariableType, AcEquationType> j,
TargetVector<AcVariableType, AcEquationType> targetVector) {
EquationSystem<AcVariableType, AcEquationType> equationSystem,
JacobianMatrix<AcVariableType, AcEquationType> j,
TargetVector<AcVariableType, AcEquationType> targetVector,
EquationVector<AcVariableType, AcEquationType> equationVector) {
this.network = Objects.requireNonNull(network);
this.parameters = Objects.requireNonNull(parameters);
this.equationSystem = Objects.requireNonNull(equationSystem);
this.j = Objects.requireNonNull(j);
this.targetVector = Objects.requireNonNull(targetVector);
this.equationVector = Objects.requireNonNull(equationVector);
}

public static List<Pair<Equation<AcVariableType, AcEquationType>, Double>> findLargestMismatches(EquationSystem<AcVariableType, AcEquationType> equationSystem, double[] mismatch, int count) {
return equationSystem.getIndex().getSortedEquationsToSolve().stream()
.map(equation -> Pair.of(equation, mismatch[equation.getColumn()]))
.filter(e -> Math.abs(e.getValue()) > Math.pow(10, -7))
.sorted(Comparator.comparingDouble((Map.Entry<Equation<AcVariableType, AcEquationType>, Double> e) -> Math.abs(e.getValue())).reversed())
.limit(count)
.collect(Collectors.toList());
}

private NewtonRaphsonStatus runIteration(double[] fx) {
private NewtonRaphsonStatus runIteration() {
LOGGER.debug("Start iteration {}", iteration);

try {
// solve f(x) = j * dx
try {
j.solveTransposed(fx);
j.solveTransposed(equationVector.getArray());
} catch (MatrixException e) {
LOGGER.error(e.toString(), e);
return NewtonRaphsonStatus.SOLVER_FAILED;
}
// f(x) now contains dx

// update x
equationSystem.getStateVector().minus(fx);
// update x and f(x) will be automatically updated
equationSystem.getStateVector().minus(equationVector.getArray());

// recalculate f(x) with new x
equationSystem.updateEquationVector(fx);

Vectors.minus(fx, targetVector.toArray());
// substract targets from f(x)
equationVector.minus(targetVector);
// f(x) now contains equation mismatches

if (LOGGER.isTraceEnabled()) {
equationSystem.findLargestMismatches(fx, 5)
findLargestMismatches(equationSystem, equationVector.getArray(), 5)
.forEach(e -> {
Equation<AcVariableType, AcEquationType> equation = e.getKey();
String elementId = equation.getElement(network).map(LfElement::getId).orElse("?");
Expand All @@ -77,7 +96,7 @@ private NewtonRaphsonStatus runIteration(double[] fx) {
}

// test stopping criteria and log norm(fx)
NewtonRaphsonStoppingCriteria.TestResult testResult = parameters.getStoppingCriteria().test(fx);
NewtonRaphsonStoppingCriteria.TestResult testResult = parameters.getStoppingCriteria().test(equationVector.getArray());

LOGGER.debug("|f(x)|={}", testResult.getNorm());

Expand Down Expand Up @@ -168,15 +187,12 @@ public NewtonRaphsonResult run(VoltageInitializer voltageInitializer) {
// initialize state vector
initStateVector(network, equationSystem, voltageInitializer);

// initialize mismatch vector (difference between equation values and targets)
double[] fx = equationSystem.createEquationVector();

Vectors.minus(fx, targetVector.toArray());
Vectors.minus(equationVector.getArray(), targetVector.getArray());

// start iterations
NewtonRaphsonStatus status = NewtonRaphsonStatus.NO_CALCULATION;
while (iteration <= parameters.getMaxIteration()) {
NewtonRaphsonStatus newStatus = runIteration(fx);
NewtonRaphsonStatus newStatus = runIteration();
if (newStatus != null) {
status = newStatus;
break;
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,7 @@
import com.powsybl.openloadflow.ac.equations.AcEquationType;
import com.powsybl.openloadflow.ac.equations.AcVariableType;
import com.powsybl.openloadflow.equations.EquationSystem;
import com.powsybl.openloadflow.equations.EquationVector;
import com.powsybl.openloadflow.equations.JacobianMatrix;
import com.powsybl.openloadflow.equations.TargetVector;
import com.powsybl.openloadflow.network.LfNetwork;
Expand All @@ -31,6 +32,8 @@ public class AcLoadFlowContext implements AutoCloseable {

private AcTargetVector targetVector;

private EquationVector<AcVariableType, AcEquationType> equationVector;

public AcLoadFlowContext(LfNetwork network, AcLoadFlowParameters parameters) {
this.network = Objects.requireNonNull(network);
this.parameters = Objects.requireNonNull(parameters);
Expand Down Expand Up @@ -65,6 +68,13 @@ public TargetVector<AcVariableType, AcEquationType> getTargetVector() {
return targetVector;
}

public EquationVector<AcVariableType, AcEquationType> getEquationVector() {
if (equationVector == null) {
equationVector = new EquationVector<>(getEquationSystem());
}
return equationVector;
}

@Override
public void close() {
if (jacobianMatrix != null) {
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -88,8 +88,12 @@ public AcLoadFlowResult run() {
voltageInitializer.prepare(context.getNetwork());

RunningContext runningContext = new RunningContext();
NewtonRaphson newtonRaphson = new NewtonRaphson(context.getNetwork(), context.getParameters().getNewtonRaphsonParameters(),
context.getEquationSystem(), context.getJacobianMatrix(), context.getTargetVector());
NewtonRaphson newtonRaphson = new NewtonRaphson(context.getNetwork(),
context.getParameters().getNewtonRaphsonParameters(),
context.getEquationSystem(),
context.getJacobianMatrix(),
context.getTargetVector(),
context.getEquationVector());

List<OuterLoop> outerLoops = context.getParameters().getOuterLoops();
List<Pair<OuterLoop, OuterLoopContextImpl>> outerLoopsAndContexts = outerLoops.stream()
Expand Down
Original file line number Diff line number Diff line change
@@ -0,0 +1,94 @@
/**
* Copyright (c) 2022, RTE (http://www.rte-france.com)
* This Source Code Form is subject to the terms of the Mozilla Public
* License, v. 2.0. If a copy of the MPL was not distributed with this
* file, You can obtain one at http://mozilla.org/MPL/2.0/.
*/
package com.powsybl.openloadflow.equations;

import java.util.Objects;

/**
* @author Geoffroy Jamgotchian <geoffroy.jamgotchian at rte-france.com>
*/
public abstract class AbstractVector<V extends Enum<V> & Quantity, E extends Enum<E> & Quantity> implements Vector {

protected final EquationSystem<V, E> equationSystem;

private double[] array;

private enum Status {
VALID,
VECTOR_INVALID, // structure has changed
VALUES_INVALID // same structure but values have to be updated
}

private Status status = Status.VECTOR_INVALID;

private final EquationSystemIndexListener<V, E> equationSystemIndexListener = new EquationSystemIndexListener<>() {

@Override
public void onEquationChange(Equation<V, E> equation, ChangeType changeType) {
invalidateVector();
}

@Override
public void onVariableChange(Variable<V> variable, ChangeType changeType) {
// nothing to do
}

@Override
public void onEquationTermChange(EquationTerm<V, E> term) {
// nothing to do
}
};

protected AbstractVector(EquationSystem<V, E> equationSystem) {
this.equationSystem = Objects.requireNonNull(equationSystem);
equationSystem.getIndex().addListener(equationSystemIndexListener);
}

protected void invalidateValues() {
if (status == Status.VALID) {
status = Status.VALUES_INVALID;
}
}

protected void invalidateVector() {
status = Status.VECTOR_INVALID;
}

protected void validate() {
status = Status.VALID;
}

@Override
public double[] getArray() {
switch (status) {
case VECTOR_INVALID:
array = createArray();
validate();
break;

case VALUES_INVALID:
updateArray(array);
validate();
break;

default:
// nothing to do
break;
}
return array;
}

protected abstract double[] createArray();

protected abstract void updateArray(double[] array);

@Override
public void minus(Vector other) {
Objects.requireNonNull(other);
Vectors.minus(getArray(), other.getArray());
}
}
Original file line number Diff line number Diff line change
Expand Up @@ -6,31 +6,23 @@
*/
package com.powsybl.openloadflow.equations;

import com.google.common.base.Stopwatch;
import com.powsybl.commons.PowsyblException;
import com.powsybl.openloadflow.network.ElementType;
import com.powsybl.openloadflow.network.LfNetwork;
import org.apache.commons.lang3.tuple.Pair;
import org.slf4j.Logger;
import org.slf4j.LoggerFactory;

import java.io.IOException;
import java.io.StringWriter;
import java.io.UncheckedIOException;
import java.io.Writer;
import java.util.*;
import java.util.concurrent.TimeUnit;
import java.util.stream.Collectors;

import static com.powsybl.openloadflow.util.Markers.PERFORMANCE_MARKER;

/**
* @author Geoffroy Jamgotchian <geoffroy.jamgotchian at rte-france.com>
*/
public class EquationSystem<V extends Enum<V> & Quantity, E extends Enum<E> & Quantity> {

private static final Logger LOGGER = LoggerFactory.getLogger(EquationSystem.class);

private final boolean indexTerms;

private final Map<Pair<Integer, E>, Equation<V, E>> equations = new HashMap<>();
Expand Down Expand Up @@ -172,27 +164,6 @@ public List<String> getColumnNames(LfNetwork network) {
.collect(Collectors.toList());
}

public double[] createEquationVector() {
double[] fx = new double[index.getSortedEquationsToSolve().size()];
updateEquationVector(fx);
return fx;
}

public void updateEquationVector(double[] fx) {
if (fx.length != index.getSortedEquationsToSolve().size()) {
throw new IllegalArgumentException("Bad equation vector length: " + fx.length);
}

Stopwatch stopwatch = Stopwatch.createStarted();

Arrays.fill(fx, 0);
for (Equation<V, E> equation : index.getSortedEquationsToSolve()) {
fx[equation.getColumn()] = equation.eval();
}

LOGGER.debug(PERFORMANCE_MARKER, "Equation vector updated in {} us", stopwatch.elapsed(TimeUnit.MICROSECONDS));
}

public void addListener(EquationSystemListener<V, E> listener) {
Objects.requireNonNull(listener);
listeners.add(listener);
Expand Down Expand Up @@ -247,13 +218,4 @@ public String writeToString(boolean writeInactiveEquations) {
public String writeToString() {
return writeToString(false);
}

public List<Pair<Equation<V, E>, Double>> findLargestMismatches(double[] mismatch, int count) {
return index.getSortedEquationsToSolve().stream()
.map(equation -> Pair.of(equation, mismatch[equation.getColumn()]))
.filter(e -> Math.abs(e.getValue()) > Math.pow(10, -7))
.sorted(Comparator.comparingDouble((Map.Entry<Equation<V, E>, Double> e) -> Math.abs(e.getValue())).reversed())
.limit(count)
.collect(Collectors.toList());
}
}
Original file line number Diff line number Diff line change
@@ -0,0 +1,60 @@
/**
* Copyright (c) 2022, RTE (http://www.rte-france.com)
* This Source Code Form is subject to the terms of the Mozilla Public
* License, v. 2.0. If a copy of the MPL was not distributed with this
* file, You can obtain one at http://mozilla.org/MPL/2.0/.
*/
package com.powsybl.openloadflow.equations;

import com.google.common.base.Stopwatch;
import org.slf4j.Logger;
import org.slf4j.LoggerFactory;

import java.util.Arrays;
import java.util.concurrent.TimeUnit;

import static com.powsybl.openloadflow.util.Markers.PERFORMANCE_MARKER;

/**
* @author Geoffroy Jamgotchian <geoffroy.jamgotchian at rte-france.com>
*/
public class EquationVector<V extends Enum<V> & Quantity, E extends Enum<E> & Quantity> extends AbstractVector<V, E>
implements StateVectorListener {

private static final Logger LOGGER = LoggerFactory.getLogger(EquationVector.class);

public EquationVector(EquationSystem<V, E> equationSystem) {
super(equationSystem);
equationSystem.getStateVector().addListener(this);
}

@Override
public void onStateUpdate() {
invalidateValues();
}

@Override
protected double[] createArray() {
double[] array = new double[equationSystem.getIndex().getSortedEquationsToSolve().size()];
updateArray(array);
return array;
}

@Override
protected void updateArray(double[] array) {
Stopwatch stopwatch = Stopwatch.createStarted();

var equations = equationSystem.getIndex().getSortedEquationsToSolve();

if (array.length != equations.size()) {
throw new IllegalArgumentException("Bad equation vector length: " + array.length);
}

Arrays.fill(array, 0); // necessary?
for (Equation<V, E> equation : equations) {
array[equation.getColumn()] = equation.eval();
}

LOGGER.debug(PERFORMANCE_MARKER, "Equation vector updated in {} us", stopwatch.elapsed(TimeUnit.MICROSECONDS));
}
}
Loading

0 comments on commit e6a160a

Please sign in to comment.