diff --git a/README.md b/README.md index 3574afb5..bad33aee 100644 --- a/README.md +++ b/README.md @@ -150,7 +150,9 @@ LoadFlow.run(network, parameters); ## Features The Knitro solver is used as a substitute for the **inner loop calculations** in the load flow process. -The **outer loops** such as distributed slack, reactive limits, etc... operate identically as when using the Newton-Raphson method. +The **outer loops** such as distributed slack, voltage monitoring, etc... operate identically as when using the Newton-Raphson method. +Note that for the reactive limits outer loop case, it is run in the standard way only when the selected Knitro solver does not directly integrate reactive +limits into the optimization model (cf [Knitro Parameters](#knitro-parameters)). ### Configuration @@ -175,6 +177,11 @@ parameters.addExtension(KnitroLoadFlowParameters.class, knitroLoadFlowParameters - **Knitro Solver Types**: - `STANDARD (default)` : the constraint satisfaction problem formulation and a direct substitute to the Newton-Raphson solver. - `RELAXED` : the optimisation problem formulation relaxing satisfaction problem. + - `USE_REACTIVE_LIMITS` : the optimization problem formulation relaxing satisfaction problem and integrating reactive limits in constraints. + Note that using this solver requires disabling the `useReactiveLimits` parameter in `LoadFlowParameters`, for compatibility reasons with the way outer + loops are launched in Open-Load-Flow. Doing so, the Open-Load-Flow checks that disable voltage control when the reactive power bounds are not sufficiently + large (see [OLF documentation](https://powsybl.readthedocs.io/projects/powsybl-open-loadflow/en/latest/loadflow/parameters.html)) are not performed, which may explain some of the differences in results between the Newton–Raphson solver and the Knitro solver. + This will be addressed in a near future. - Use `setKnitroSolverType` in the `KnitroLoadFlowParameters` extension. 2. **Voltage Bounds**: diff --git a/src/main/java/com/powsybl/openloadflow/knitro/solver/AbstractKnitroProblem.java b/src/main/java/com/powsybl/openloadflow/knitro/solver/AbstractKnitroProblem.java index 3b3950b7..ca3c212d 100644 --- a/src/main/java/com/powsybl/openloadflow/knitro/solver/AbstractKnitroProblem.java +++ b/src/main/java/com/powsybl/openloadflow/knitro/solver/AbstractKnitroProblem.java @@ -31,7 +31,7 @@ * - Initialization and definition of variable bounds for the optimization problem. * - Definition of constraints (including those evaluated via callbacks in {@link KnitroCallbacks}). * - Representation of the constraint Jacobian for the problem. - * This class can be extended to customize any of these features (e.g., in {@link RelaxedKnitroSolver.RelaxedKnitroProblem}). + * This class can be extended to customize any of these features (e.g., in {@link AbstractRelaxedKnitroSolver.AbstractRelaxedKnitroProblem}). * For example, if you modify the optimization problem, you may also need to update the initialization of additional variables. * * @author Jeanne Archambault {@literal } @@ -51,17 +51,26 @@ public abstract class AbstractKnitroProblem extends KNProblem { protected final int numberOfPowerFlowVariables; protected List> activeConstraints = new ArrayList<>(); protected final List nonlinearConstraintIndexes = new ArrayList<>(); + protected final int numTotalVariables; protected AbstractKnitroProblem(LfNetwork network, EquationSystem equationSystem, TargetVector targetVector, JacobianMatrix jacobianMatrix, - KnitroSolverParameters knitroParameters, int numTotalVariables) { - super(numTotalVariables, equationSystem.getIndex().getSortedEquationsToSolve().size()); + KnitroSolverParameters knitroParameters) { + this(network, equationSystem, targetVector, jacobianMatrix, knitroParameters, 0, 0); + } + + protected AbstractKnitroProblem(LfNetwork network, EquationSystem equationSystem, + TargetVector targetVector, JacobianMatrix jacobianMatrix, + KnitroSolverParameters knitroParameters, int numAdditionalVariables, int numAdditionalConstraints) { + super(equationSystem.getIndex().getSortedVariablesToFind().size() + numAdditionalVariables, + equationSystem.getIndex().getSortedEquationsToSolve().size() + numAdditionalConstraints); + this.numberOfPowerFlowVariables = equationSystem.getIndex().getSortedVariablesToFind().size(); + this.numTotalVariables = numberOfPowerFlowVariables + numAdditionalVariables; this.network = network; this.equationSystem = equationSystem; this.targetVector = targetVector; this.jacobianMatrix = jacobianMatrix; this.knitroParameters = knitroParameters; - this.numberOfPowerFlowVariables = equationSystem.getIndex().getSortedVariablesToFind().size(); } /** @@ -69,9 +78,8 @@ protected AbstractKnitroProblem(LfNetwork network, EquationSystem variableTypes = new ArrayList<>(Collections.nCopies(numTotalVariables, KNConstants.KN_VARTYPE_CONTINUOUS)); List lowerBounds = new ArrayList<>(Collections.nCopies(numTotalVariables, -KNConstants.KN_INFINITY)); List upperBounds = new ArrayList<>(Collections.nCopies(numTotalVariables, KNConstants.KN_INFINITY)); @@ -95,7 +103,10 @@ protected void initializeVariables(VoltageInitializer voltageInitializer, int nu } // Allow subclasses to modify bounds and initial values (e.g., for slack variables) - initializeCustomizedVariables(lowerBounds, upperBounds, initialValues, numTotalVariables); + initializeCustomizedVariables(lowerBounds, upperBounds, initialValues); + + // Set up scaling factors + setUpScalingFactors(); setVarLoBnds(lowerBounds); setVarUpBnds(upperBounds); @@ -110,10 +121,18 @@ protected void initializeVariables(VoltageInitializer voltageInitializer, int nu * @param lowerBounds Lower bounds list to modify. * @param upperBounds Upper bounds list to modify. * @param initialValues Initial values list to modify. - * @param numTotalVariables Total number of variables. */ protected void initializeCustomizedVariables(List lowerBounds, List upperBounds, - List initialValues, int numTotalVariables) { + List initialValues) { + // no customization by default + } + + /** + * Allows subclasses to utilize scaling. + * Default implementation does nothing. + * + */ + protected void setUpScalingFactors() throws KNException { // no customization by default } @@ -131,7 +150,7 @@ protected void setupConstraints() throws KNException { NonLinearExternalSolverUtils solverUtils = new NonLinearExternalSolverUtils(); // add linear constraints and fill the list of non-linear constraints - addLinearConstraints(activeConstraints, solverUtils, nonlinearConstraintIndexes); + addLinearConstraints(activeConstraints, solverUtils); // pass to Knitro the indexes of non-linear constraints, that will be evaluated in the callback function setMainCallbackCstIndexes(nonlinearConstraintIndexes); @@ -145,14 +164,12 @@ protected void setupConstraints() throws KNException { * * @param sortedEquationsToSolve Sorted list of equations to solve. * @param solverUtils Utilities to extract linear constraints. - * @param nonLinearConstraintIds Output list of indices of non-linear constraints. */ protected void addLinearConstraints(List> sortedEquationsToSolve, - NonLinearExternalSolverUtils solverUtils, - List nonLinearConstraintIds) { + NonLinearExternalSolverUtils solverUtils) { for (int equationId = 0; equationId < sortedEquationsToSolve.size(); equationId++) { - addConstraint(equationId, sortedEquationsToSolve, solverUtils, nonLinearConstraintIds); + addConstraint(equationId, sortedEquationsToSolve, solverUtils); } } @@ -163,10 +180,9 @@ protected void addLinearConstraints(List> sortedEquationsToSolve, - NonLinearExternalSolverUtils solverUtils, List nonLinearConstraintIds) { + NonLinearExternalSolverUtils solverUtils) { Equation equation = sortedEquationsToSolve.get(equationId); AcEquationType equationType = equation.getType(); @@ -190,7 +206,7 @@ protected void addConstraint(int equationId, List} diff --git a/src/main/java/com/powsybl/openloadflow/knitro/solver/AbstractRelaxedKnitroSolver.java b/src/main/java/com/powsybl/openloadflow/knitro/solver/AbstractRelaxedKnitroSolver.java new file mode 100644 index 00000000..55f4e1ac --- /dev/null +++ b/src/main/java/com/powsybl/openloadflow/knitro/solver/AbstractRelaxedKnitroSolver.java @@ -0,0 +1,446 @@ +/** + * Copyright (c) 2025, Artelys (http://www.artelys.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/. + * SPDX-License-Identifier: MPL-2.0 + */ +package com.powsybl.openloadflow.knitro.solver; + +import com.artelys.knitro.api.*; +import com.artelys.knitro.api.callbacks.KNEvalGACallback; +import com.powsybl.commons.PowsyblException; +import com.powsybl.openloadflow.ac.equations.AcEquationType; +import com.powsybl.openloadflow.ac.equations.AcVariableType; +import com.powsybl.openloadflow.equations.*; +import com.powsybl.openloadflow.network.LfBus; +import com.powsybl.openloadflow.network.LfNetwork; +import com.powsybl.openloadflow.util.PerUnit; +import org.slf4j.Logger; +import org.slf4j.LoggerFactory; + +import java.util.*; + +/** + * Abstract class for relaxed Knitro solvers, solving the open load flow equation system by minimizing constraint violations through relaxation. + * It provides common functionality, including: + * - Post-processing of the computed solutions, including the reporting of relaxation variables. + * - An extended optimization problem formulation dedicated to solving the open load flow equation system, including relaxations. + * This class can be extended to add custom behavior to any of these features (e.g., in {@link com.powsybl.openloadflow.knitro.solver.UseReactiveLimitsKnitroSolver}). + * For example, if you modify the optimization problem, you may also need to update the solution-processing logic. + * + * @author Martin Debouté {@literal } + * @author Amine Makhen {@literal } + * @author Pierre Arvy {@literal } + */ +public abstract class AbstractRelaxedKnitroSolver extends AbstractKnitroSolver { + + private static final Logger LOGGER = LoggerFactory.getLogger(AbstractRelaxedKnitroSolver.class); + + // Penalty weights in the objective function + protected static final double WEIGHT_P_PENAL = 1.0; + protected static final double WEIGHT_Q_PENAL = 1.0; + protected static final double WEIGHT_V_PENAL = 1.0; + + // Weights of the linear in the objective function + protected static final double WEIGHT_ABSOLUTE_PENAL = 3.0; + + // Total number of variables (including power flow and slack variables) + protected int numSlackVariables; + + // Number of equations for active power (P), reactive power (Q), and voltage magnitude (V) + protected final int numPEquations; + protected final int numQEquations; + protected final int numVEquations; + + // Starting indices for slack variables in the variable vector + protected final int slackPStartIndex; + protected final int slackQStartIndex; + protected final int slackVStartIndex; + + // Mappings from global equation indices to local indices by equation type + protected final Map pEquationLocalIds; + protected final Map qEquationLocalIds; + protected final Map vEquationLocalIds; + + protected AbstractRelaxedKnitroSolver(LfNetwork network, KnitroSolverParameters knitroParameters, EquationSystem equationSystem, + JacobianMatrix j, TargetVector targetVector, + EquationVector equationVector, boolean detailedReport) { + super(network, knitroParameters, equationSystem, j, targetVector, equationVector, detailedReport); + + List> sortedEquations = equationSystem.getIndex().getSortedEquationsToSolve(); + + // Count number of equations by type + this.numPEquations = (int) sortedEquations.stream().filter(e -> e.getType() == AcEquationType.BUS_TARGET_P).count(); + this.numQEquations = (int) sortedEquations.stream().filter(e -> e.getType() == AcEquationType.BUS_TARGET_Q).count(); + this.numVEquations = (int) sortedEquations.stream().filter(e -> e.getType() == AcEquationType.BUS_TARGET_V).count(); + + this.numSlackVariables = 2 * (numPEquations + numQEquations + numVEquations); + + // the slack variables start after power flow variables + this.slackPStartIndex = equationSystem.getIndex().getSortedVariablesToFind().size(); + this.slackQStartIndex = slackPStartIndex + 2 * numPEquations; + this.slackVStartIndex = slackQStartIndex + 2 * numQEquations; + + // Map equations to local indices + this.pEquationLocalIds = new HashMap<>(); + this.qEquationLocalIds = new HashMap<>(); + this.vEquationLocalIds = new HashMap<>(); + + int pCounter = 0; + int qCounter = 0; + int vCounter = 0; + + for (int i = 0; i < sortedEquations.size(); i++) { + AcEquationType type = sortedEquations.get(i).getType(); + + switch (type) { + case BUS_TARGET_P -> pEquationLocalIds.put(i, pCounter++); + case BUS_TARGET_Q -> qEquationLocalIds.put(i, qCounter++); + case BUS_TARGET_V -> vEquationLocalIds.put(i, vCounter++); + default -> { + // Other equation types don't require slack variables + } + } + } + } + + @Override + protected void processSolution(KNSolver solver, KNSolution solution, KNProblem problemInstance) { + super.processSolution(solver, solution, problemInstance); + + List x = solution.getX(); + + // ========== Slack Logging ========== + logSlackValues("P", slackPStartIndex, numPEquations, x); + logSlackValues("Q", slackQStartIndex, numQEquations, x); + logSlackValues("V", slackVStartIndex, numVEquations, x); + + // ========== Penalty Computation ========== + double penaltyP = computeSlackPenalty(x, slackPStartIndex, numPEquations, WEIGHT_P_PENAL); + double penaltyQ = computeSlackPenalty(x, slackQStartIndex, numQEquations, WEIGHT_Q_PENAL); + double penaltyV = computeSlackPenalty(x, slackVStartIndex, numVEquations, WEIGHT_V_PENAL); + double totalPenalty = penaltyP + penaltyQ + penaltyV; + + LOGGER.info("==== Slack penalty details ===="); + LOGGER.info("Penalty P = {}", penaltyP); + LOGGER.info("Penalty Q = {}", penaltyQ); + LOGGER.info("Penalty V = {}", penaltyV); + LOGGER.info("Total penalty = {}", totalPenalty); + } + + /** + * Logs information like bus name and slack value for most significant slack variables. + * + * @param type The slack variable type. + * @param startIndex The start index of slack variables associated to the given type. + * @param count The maximum number of slack variables associated to the given type. + * @param x The variable values as returned by solver. + */ + protected void logSlackValues(String type, int startIndex, int count, List x) { + LOGGER.debug("==== Slack diagnostics for {} (p.u. and physical units) ====", type); + + for (int i = 0; i < count; i++) { + double sm = x.get(startIndex + 2 * i); + double sp = x.get(startIndex + 2 * i + 1); + double epsilon = sp - sm; + + // Get significant slack values above threshold + boolean shouldSkip = Math.abs(epsilon) <= knitroParameters.getSlackThreshold(); + String name = null; + String interpretation = null; + + if (!shouldSkip) { + name = getSlackVariableBusName(i, type); + + switch (type) { + case "P" -> interpretation = String.format("ΔP = %.4f p.u. (%.1f MW)", epsilon, epsilon * PerUnit.SB); + case "Q" -> interpretation = String.format("ΔQ = %.4f p.u. (%.1f MVAr)", epsilon, epsilon * PerUnit.SB); + case "V" -> { + var bus = network.getBusById(name); + if (bus == null) { + LOGGER.warn("Bus {} not found while logging V slack.", name); + shouldSkip = true; + } else { + interpretation = String.format("ΔV = %.4f p.u. (%.1f kV)", epsilon, epsilon * bus.getNominalV()); + } + } + default -> interpretation = "Unknown slack type"; + } + } + + if (shouldSkip) { + continue; + } + + String msg = String.format("Slack %s[ %s ] → Sm = %.4f, Sp = %.4f → %s", type, name, sm, sp, interpretation); + LOGGER.debug(msg); + } + } + + /** + * Finds the bus associated to a slack variable. + * + * @param index The index of the slack variable + * @param type The slack variable type. + * @return The id of the bus associated to the slack variable. + */ + private String getSlackVariableBusName(Integer index, String type) { + Set> equationSet = switch (type) { + case "P" -> pEquationLocalIds.entrySet(); + case "Q" -> qEquationLocalIds.entrySet(); + case "V" -> vEquationLocalIds.entrySet(); + default -> throw new IllegalStateException("Unexpected variable type: " + type); + }; + + Optional varIndexOptional = equationSet.stream() + .filter(entry -> index.equals(entry.getValue())) + .map(Map.Entry::getKey) + .findAny(); + + int varIndex; + if (varIndexOptional.isPresent()) { + varIndex = varIndexOptional.get(); + } else { + throw new PowsyblException("Variable index associated with slack variable " + type + " was not found"); + } + + LfBus bus = network.getBus(equationSystem.getIndex().getSortedEquationsToSolve().get(varIndex).getElementNum()); + + return bus.getId(); + } + + /** + * Calculates the total loss associated to a slack variable type + * + * @param x The variable values as returned by solver. + * @param startIndex The start index of slack variables associated to the given type. + * @param count The maximum number of slack variables associated to the given type. + * @param weight The weight inf front of the given slack variables terms + * @return The total penalty associated to the slack variables type. + */ + double computeSlackPenalty(List x, int startIndex, int count, double weight) { + double penalty = 0.0; + for (int i = 0; i < count; i++) { + double sm = x.get(startIndex + 2 * i); + double sp = x.get(startIndex + 2 * i + 1); + double diff = sp - sm; + penalty += weight * (diff * diff); // Quadratic terms + penalty += weight * WEIGHT_ABSOLUTE_PENAL * (sp + sm); // Linear terms + } + return penalty; + } + + /** + * Optimization problem-solving the open load flow equation system by minimizing constraint violations through relaxation. + */ + public abstract class AbstractRelaxedKnitroProblem extends AbstractKnitroProblem { + + /** + * Relaxed Knitro problem definition including: + * - initialization of variables (types, bounds, initial state) + * - definition of linear constraints + * - definition of non-linear constraints, evaluated in extended the callback function + * - definition of the extended Jacobian matrix passed to Knitro to solve the problem + * - definition of the objective function to be minimized (equation system violations) + */ + protected AbstractRelaxedKnitroProblem(LfNetwork network, EquationSystem equationSystem, + TargetVector targetVector, JacobianMatrix jacobianMatrix, + KnitroSolverParameters knitroParameters, int numAdditionalVariables, int numAdditionalConstraints) { + super(network, equationSystem, targetVector, jacobianMatrix, knitroParameters, numAdditionalVariables, numAdditionalConstraints); + } + + void addObjectiveFunction(int numPEquations, int slackPStartIndex, int numQEquations, int slackQStartIndex, + int numVEquations, int slackVStartIndex) throws KNException { + // initialise lists to track quadratic objective function terms of the form: a * x1 * x2 + List quadRows = new ArrayList<>(); // list of indexes of the first variable x1 + List quadCols = new ArrayList<>(); // list of indexes of the second variable x2 + List quadCoefs = new ArrayList<>(); // list of indexes of the coefficient a + + // initialise lists to track linear objective function terms of the form: a * x + List linIndexes = new ArrayList<>(); // list of indexes of the variable x + List linCoefs = new ArrayList<>(); // list of indexes of the coefficient a + + // add slack penalty terms, for each slack type, of the form: (Sp - Sm)^2 = Sp^2 + Sm^2 - 2*Sp*Sm + linear terms from the absolute value + addSlackObjectiveTerms(numPEquations, slackPStartIndex, AbstractRelaxedKnitroSolver.WEIGHT_P_PENAL, AbstractRelaxedKnitroSolver.WEIGHT_ABSOLUTE_PENAL, quadRows, quadCols, quadCoefs, linIndexes, linCoefs); + addSlackObjectiveTerms(numQEquations, slackQStartIndex, AbstractRelaxedKnitroSolver.WEIGHT_Q_PENAL, AbstractRelaxedKnitroSolver.WEIGHT_ABSOLUTE_PENAL, quadRows, quadCols, quadCoefs, linIndexes, linCoefs); + addSlackObjectiveTerms(numVEquations, slackVStartIndex, AbstractRelaxedKnitroSolver.WEIGHT_V_PENAL, AbstractRelaxedKnitroSolver.WEIGHT_ABSOLUTE_PENAL, quadRows, quadCols, quadCoefs, linIndexes, linCoefs); + + setObjectiveQuadraticPart(quadRows, quadCols, quadCoefs); + setObjectiveLinearPart(linIndexes, linCoefs); + } + + /** + * Adds quadratic and linear terms related to slack variables to the objective function. + */ + void addSlackObjectiveTerms(int numEquations, int slackStartIdx, double weight, double lambda, + List quadRows, List quadCols, List quadCoefs, + List linIndexes, List linCoefs) { + for (int i = 0; i < numEquations; i++) { + int idxSm = slackStartIdx + 2 * i; // negative slack variable index + int idxSp = slackStartIdx + 2 * i + 1; // positive slack variable index + + // Add quadratic terms: weight * (sp^2 + sm^2 - 2 * sp * sm) + + // add first quadratic term : weight * sp^2 + quadRows.add(idxSp); + quadCols.add(idxSp); + quadCoefs.add(weight); + + // add second quadratic term : weight * sm^2 + quadRows.add(idxSm); + quadCols.add(idxSm); + quadCoefs.add(weight); + + // add third quadratic term : weight * (- 2 * sp * sm) + quadRows.add(idxSp); + quadCols.add(idxSm); + quadCoefs.add(-2 * weight); + + // Add linear terms: weight * lambda * (sp + sm) + + // add first linear term : weight * lambda * sp + linIndexes.add(idxSp); + linCoefs.add(lambda * weight); + + // add second linear term : weight * lambda * sm + linIndexes.add(idxSm); + linCoefs.add(lambda * weight); + } + } + + @Override + protected void initializeCustomizedVariables(List lowerBounds, List upperBounds, + List initialValues) { + // set a lower bound to slack variables (>= 0) + // initial values have already been set to 0 + for (int i = numberOfPowerFlowVariables; i < numTotalVariables; i++) { + lowerBounds.set(i, 0.0); + } + } + + @Override + protected void addAdditionalConstraintVariables(int equationId, AcEquationType equationType, + List varIndices, List coefficients) { + // Add slack variables if applicable + int slackBase = getSlackIndexBase(equationType, equationId); + if (slackBase >= 0) { + varIndices.add(slackBase); // Sm + varIndices.add(slackBase + 1); // Sp + coefficients.add(1.0); + coefficients.add(-1.0); + } + } + + @Override + protected void addAdditionalJacobianVariables(int constraintIndex, + Equation equation, + List variableIndices) { + AcEquationType equationType = equation.getType(); + // get slack variable local index (within its equation type) + int slackStart = switch (equationType) { + case BUS_TARGET_P -> pEquationLocalIds.getOrDefault(constraintIndex, -1); + case BUS_TARGET_Q -> qEquationLocalIds.getOrDefault(constraintIndex, -1); + case BUS_TARGET_V -> vEquationLocalIds.getOrDefault(constraintIndex, -1); + default -> -1; + }; + + if (slackStart >= 0) { + // get slack variable type starting index (within total variables' indexes) + int slackBaseIndex = switch (equationType) { + case BUS_TARGET_P -> slackPStartIndex; + case BUS_TARGET_Q -> slackQStartIndex; + case BUS_TARGET_V -> slackVStartIndex; + default -> throw new IllegalStateException("Unexpected constraint type: " + equationType); + }; + // get slack variables Sm and Sp indexes + variableIndices.add(slackBaseIndex + 2 * slackStart); // Sm + variableIndices.add(slackBaseIndex + 2 * slackStart + 1); // Sp + } + } + + /** + * Returns the base index of the slack variable associated with a given equation type and ID. + * + * @param equationType Type of the equation (P, Q, or V). + * @param equationId Index of the equation. + * @return Base index of the corresponding slack variable, or -1 if not applicable. + */ + protected int getSlackIndexBase(AcEquationType equationType, int equationId) { + return switch (equationType) { + case BUS_TARGET_P -> pEquationLocalIds.getOrDefault(equationId, -1) >= 0 + ? slackPStartIndex + 2 * pEquationLocalIds.get(equationId) : -1; + case BUS_TARGET_Q -> qEquationLocalIds.getOrDefault(equationId, -1) >= 0 + ? slackQStartIndex + 2 * qEquationLocalIds.get(equationId) : -1; + case BUS_TARGET_V -> vEquationLocalIds.getOrDefault(equationId, -1) >= 0 + ? slackVStartIndex + 2 * vEquationLocalIds.get(equationId) : -1; + default -> -1; + }; + } + + @Override + protected KNEvalGACallback createGradientCallback(JacobianMatrix jacobianMatrix, + List listNonZerosCtsDense, List listNonZerosVarsDense, + List listNonZerosCtsSparse, List listNonZerosVarsSparse) { + + return new AbstractRelaxedKnitroProblem.RelaxedCallbackEvalG(jacobianMatrix, listNonZerosCtsDense, listNonZerosVarsDense, + listNonZerosCtsSparse, listNonZerosVarsSparse, network, + equationSystem, knitroParameters, numberOfPowerFlowVariables); + } + + /** + * Callback used by Knitro to evaluate the non-linear parts of the objective and constraint functions. + */ + public static class RelaxedCallbackEvalFC extends KnitroCallbacks.BaseCallbackEvalFC { + + private final AbstractRelaxedKnitroProblem problemInstance; + + RelaxedCallbackEvalFC(AbstractRelaxedKnitroProblem problemInstance, + List> sortedEquationsToSolve, + List nonLinearConstraintIds) { + super(sortedEquationsToSolve, nonLinearConstraintIds); + this.problemInstance = problemInstance; + } + + @Override + protected double addModificationOfNonLinearConstraints(int equationId, AcEquationType equationType, + List x) { + int slackIndexBase = problemInstance.getSlackIndexBase(equationType, equationId); + if (slackIndexBase >= 0) { + double sm = x.get(slackIndexBase); // negative slack + double sp = x.get(slackIndexBase + 1); // positive slack + return sp - sm; // add slack contribution + } + return 0; + } + } + + /** + * Callback used by Knitro to evaluate the gradient (Jacobian matrix) of the constraints. + * Only constraints (no objective) are handled here. + */ + public static class RelaxedCallbackEvalG extends KnitroCallbacks.BaseCallbackEvalG { + + RelaxedCallbackEvalG(JacobianMatrix jacobianMatrix, + List denseConstraintIndices, List denseVariableIndices, + List sparseConstraintIndices, List sparseVariableIndices, + LfNetwork network, EquationSystem equationSystem, + KnitroSolverParameters knitroParameters, int numLFVariables) { + + super(jacobianMatrix, denseConstraintIndices, denseVariableIndices, sparseConstraintIndices, sparseVariableIndices, + network, equationSystem, knitroParameters, numLFVariables); + } + + @Override + protected double computeModifiedJacobianValue(int variableIndex, int constraintIndex) { + if (((variableIndex - numLFVariables) & 1) == 0) { + // set Jacobian entry to -1.0 if slack variable is Sm + return -1.0; + } else { + // set Jacobian entry to +1.0 if slack variable is Sp + return 1.0; + } + } + } + } +} diff --git a/src/main/java/com/powsybl/openloadflow/knitro/solver/KnitroCallbacks.java b/src/main/java/com/powsybl/openloadflow/knitro/solver/KnitroCallbacks.java index 3d4c698c..799724e1 100644 --- a/src/main/java/com/powsybl/openloadflow/knitro/solver/KnitroCallbacks.java +++ b/src/main/java/com/powsybl/openloadflow/knitro/solver/KnitroCallbacks.java @@ -19,6 +19,7 @@ import org.slf4j.Logger; import org.slf4j.LoggerFactory; +import java.util.Arrays; import java.util.List; import static com.google.common.primitives.Doubles.toArray; @@ -113,7 +114,7 @@ public void evaluateFC(final List x, final List obj, final List< * Allows modification of the constraint value. * Default implementation returns no modification. * Should be overridden by subclasses that modify the callbacks of Jacobian matrix - * (e.g., to add relaxation variables in {@link RelaxedKnitroSolver.RelaxedKnitroProblem}). + * (e.g., to add relaxation variables in {@link AbstractRelaxedKnitroSolver.AbstractRelaxedKnitroProblem}). * * @param equationId The equation ID. * @param equationType The equation type. @@ -212,27 +213,30 @@ protected void fillJacobianValues(List constraintIndices, List try { double value; - // Find matching (variableIndex, constraintIndex) entry in sparse column - int colStart = columnStart[constraintIndex]; - - // Check if we moved to the next constraint - if (currentConstraint != constraintIndex) { - // In this case, restart variable tracking slacks indices to 0 and change the currently handled constraint - iRowIndices = 0; - currentConstraint = constraintIndex; - } - - // Check if current variableIndex is a slack variable index - if (variableIndex >= numLFVariables) { - // In this case, add corresponding value - value = computeModifiedJacobianValue(variableIndex, constraintIndex); - // When in dense mode, the constructed index list and the index list from powsybl don't match - } else if (rowIndices[colStart + iRowIndices] != variableIndex) { - // In this case, set corresponding value to zero - value = 0.0; + if (constraintIndex < Arrays.stream(columnStart).count()) { // Find matching (variableIndex, constraintIndex) entry in sparse column + int colStart = columnStart[constraintIndex]; + + // Check if we moved to the next constraint + if (currentConstraint != constraintIndex) { + // In this case, restart variable tracking slacks indices to 0 and change the currently handled constraint + iRowIndices = 0; + currentConstraint = constraintIndex; + } + + // Check if current variableIndex is a slack variable index + if (variableIndex >= numLFVariables) { + // In this case, add corresponding value + value = computeModifiedJacobianValue(variableIndex, constraintIndex); + // When in dense mode, the constructed index list and the index list from powsybl don't match + } else if (rowIndices[colStart + iRowIndices] != variableIndex) { + // In this case, set corresponding value to zero + value = 0.0; + } else { + // This is the case of a LF variable with a non-zero coefficient in the Jacobian + value = values[colStart + iRowIndices++]; + } } else { - // This is the case of a LF variable with a non-zero coefficient in the Jacobian - value = values[colStart + iRowIndices++]; + value = computeAddedConstraintJacobianValue(variableIndex, constraintIndex); } jac.set(index, value); @@ -247,11 +251,20 @@ protected void fillJacobianValues(List constraintIndices, List /** * Computes a modified Jacobian value. Default implementation returns 0. * Should be overridden by subclasses that modify the callbacks of Jacobian matrix - * (e.g., to add relaxation variables in {@link RelaxedKnitroSolver.RelaxedKnitroProblem}). + * (e.g., to add relaxation variables in {@link AbstractRelaxedKnitroSolver.AbstractRelaxedKnitroProblem}). */ protected double computeModifiedJacobianValue(int variableIndex, int constraintIndex) { return 0.0; } + + /** + * Computes the value of an added constraint. Default implementation returns 0. + * Should be overridden by subclasses that modify the callbacks of Jacobian matrix + * (e.g., to add relaxation variables in {@link UseReactiveLimitsKnitroSolver.UseReactiveLimitsKnitroProblem}). + */ + protected double computeAddedConstraintJacobianValue(int variableIndex, int constraintIndex) { + return 0.0; + } } } diff --git a/src/main/java/com/powsybl/openloadflow/knitro/solver/KnitroSolver.java b/src/main/java/com/powsybl/openloadflow/knitro/solver/KnitroSolver.java index d7303069..9d1f285d 100644 --- a/src/main/java/com/powsybl/openloadflow/knitro/solver/KnitroSolver.java +++ b/src/main/java/com/powsybl/openloadflow/knitro/solver/KnitroSolver.java @@ -70,13 +70,12 @@ private KnitroProblem(LfNetwork lfNetwork, EquationSystem new KnitroSolver(network, knitroSolverParameters, equationSystem, j, targetVector, equationVector, parameters.isDetailedReport()); case RELAXED -> new RelaxedKnitroSolver(network, knitroSolverParameters, equationSystem, j, targetVector, equationVector, parameters.isDetailedReport()); + case USE_REACTIVE_LIMITS -> new UseReactiveLimitsKnitroSolver(network, knitroSolverParameters, equationSystem, j, targetVector, equationVector, parameters.isDetailedReport()); }; } @Override public void checkSolverAndParameterConsistency(LoadFlowParameters loadFlowParameters, OpenLoadFlowParameters openLoadFlowParameters) { - // no current incompatibilities between Knitro Solver and parameters + KnitroLoadFlowParameters knitroLoadFlowParameters = loadFlowParameters.getExtension(KnitroLoadFlowParameters.class); + if (knitroLoadFlowParameters != null && knitroLoadFlowParameters.getKnitroSolverType() == KnitroSolverParameters.SolverType.USE_REACTIVE_LIMITS) { + // since reactive limits are taken into account in the Knitro solver, there is no need to activate the outer loop + if (loadFlowParameters.isUseReactiveLimits()) { + throw new PowsyblException("Knitro generator reactive limits and reactive limits outer loop cannot work simultaneously: useReactiveLimits LoadFlowParameter should be switched to false"); + } + + // exact dense Jacobian computation mode is not supported by Knitro generator reactive limits solver + if (knitroLoadFlowParameters.getGradientComputationMode() == 1 && knitroLoadFlowParameters.getGradientUserRoutine() == 1) { + throw new PowsyblException("Knitro generator reactive limits is incompatible with exact dense jacobian computation mode: gradientUserRoutine KnitroLoadFlowParameters should be switched to 1"); + } + } } } diff --git a/src/main/java/com/powsybl/openloadflow/knitro/solver/KnitroSolverParameters.java b/src/main/java/com/powsybl/openloadflow/knitro/solver/KnitroSolverParameters.java index 862dd415..c7ad4f46 100644 --- a/src/main/java/com/powsybl/openloadflow/knitro/solver/KnitroSolverParameters.java +++ b/src/main/java/com/powsybl/openloadflow/knitro/solver/KnitroSolverParameters.java @@ -295,6 +295,7 @@ public String toString() { public enum SolverType { STANDARD, - RELAXED + RELAXED, + USE_REACTIVE_LIMITS, } } diff --git a/src/main/java/com/powsybl/openloadflow/knitro/solver/RelaxedKnitroSolver.java b/src/main/java/com/powsybl/openloadflow/knitro/solver/RelaxedKnitroSolver.java index e978ee94..58618905 100644 --- a/src/main/java/com/powsybl/openloadflow/knitro/solver/RelaxedKnitroSolver.java +++ b/src/main/java/com/powsybl/openloadflow/knitro/solver/RelaxedKnitroSolver.java @@ -7,58 +7,25 @@ */ package com.powsybl.openloadflow.knitro.solver; -import com.artelys.knitro.api.*; -import com.artelys.knitro.api.callbacks.KNEvalGACallback; +import com.artelys.knitro.api.KNException; +import com.artelys.knitro.api.KNProblem; import com.powsybl.commons.PowsyblException; import com.powsybl.openloadflow.ac.equations.AcEquationType; import com.powsybl.openloadflow.ac.equations.AcVariableType; -import com.powsybl.openloadflow.equations.*; -import com.powsybl.openloadflow.network.LfBus; +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; import com.powsybl.openloadflow.network.util.VoltageInitializer; -import org.slf4j.Logger; -import org.slf4j.LoggerFactory; - -import java.util.*; /** * Relaxed Knitro solver, solving the open load flow equation system by minimizing constraint violations through relaxation. - * This class additionally provides: - * - Post-processing of the computed solutions, including the reporting of relaxation variables. - * - An extended optimization problem formulation dedicated to solving the open load flow equation system. * * @author Martin Debouté {@literal } * @author Amine Makhen {@literal } */ -public class RelaxedKnitroSolver extends AbstractKnitroSolver { - - private static final Logger LOGGER = LoggerFactory.getLogger(RelaxedKnitroSolver.class); - - // Penalty weights in the objective function - private static final double WEIGHT_P_PENAL = 1.0; - private static final double WEIGHT_Q_PENAL = 1.0; - private static final double WEIGHT_V_PENAL = 1.0; - - // Weights of the linear in the objective function - private static final double WEIGHT_ABSOLUTE_PENAL = 3.0; - - // Total number of variables (including power flow and slack variables) - private final int numberOfVariables; - - // Number of equations for active power (P), reactive power (Q), and voltage magnitude (V) - private final int numPEquations; - private final int numQEquations; - private final int numVEquations; - - // Starting indices for slack variables in the variable vector - private final int slackPStartIndex; - private final int slackQStartIndex; - private final int slackVStartIndex; - - // Mappings from global equation indices to local indices by equation type - private final Map pEquationLocalIds; - private final Map qEquationLocalIds; - private final Map vEquationLocalIds; +public class RelaxedKnitroSolver extends AbstractRelaxedKnitroSolver { public RelaxedKnitroSolver( LfNetwork network, @@ -70,45 +37,6 @@ public RelaxedKnitroSolver( boolean detailedReport) { super(network, knitroParameters, equationSystem, j, targetVector, equationVector, detailedReport); - - // Number of variables in the equations system of open load flow - int numberOfPowerFlowVariables = equationSystem.getIndex().getSortedVariablesToFind().size(); - - List> sortedEquations = equationSystem.getIndex().getSortedEquationsToSolve(); - - // Count number of equations by type - this.numPEquations = (int) sortedEquations.stream().filter(e -> e.getType() == AcEquationType.BUS_TARGET_P).count(); - this.numQEquations = (int) sortedEquations.stream().filter(e -> e.getType() == AcEquationType.BUS_TARGET_Q).count(); - this.numVEquations = (int) sortedEquations.stream().filter(e -> e.getType() == AcEquationType.BUS_TARGET_V).count(); - - int numSlackVariables = 2 * (numPEquations + numQEquations + numVEquations); - this.numberOfVariables = numberOfPowerFlowVariables + numSlackVariables; - - this.slackPStartIndex = numberOfPowerFlowVariables; - this.slackQStartIndex = slackPStartIndex + 2 * numPEquations; - this.slackVStartIndex = slackQStartIndex + 2 * numQEquations; - - // Map equations to local indices - this.pEquationLocalIds = new HashMap<>(); - this.qEquationLocalIds = new HashMap<>(); - this.vEquationLocalIds = new HashMap<>(); - - int pCounter = 0; - int qCounter = 0; - int vCounter = 0; - - for (int i = 0; i < sortedEquations.size(); i++) { - AcEquationType type = sortedEquations.get(i).getType(); - - switch (type) { - case BUS_TARGET_P -> pEquationLocalIds.put(i, pCounter++); - case BUS_TARGET_Q -> qEquationLocalIds.put(i, qCounter++); - case BUS_TARGET_V -> vEquationLocalIds.put(i, vCounter++); - default -> { - // Other equation types don't require slack variables - } - } - } } @Override @@ -119,167 +47,24 @@ public String getName() { @Override protected KNProblem createKnitroProblem(VoltageInitializer voltageInitializer) { try { - return new RelaxedKnitroProblem(network, equationSystem, targetVector, j, voltageInitializer, knitroParameters); + return new RelaxedKnitroProblem(network, equationSystem, targetVector, j, knitroParameters, voltageInitializer); } catch (KNException e) { throw new PowsyblException("Failed to create relaxed Knitro problem", e); } } - @Override - protected void processSolution(KNSolver solver, KNSolution solution, KNProblem problemInstance) { - super.processSolution(solver, solution, problemInstance); - - List x = solution.getX(); - - // ========== Slack Logging ========== - logSlackValues("P", slackPStartIndex, numPEquations, x); - logSlackValues("Q", slackQStartIndex, numQEquations, x); - logSlackValues("V", slackVStartIndex, numVEquations, x); - - // ========== Penalty Computation ========== - double penaltyP = computeSlackPenalty(x, slackPStartIndex, numPEquations, WEIGHT_P_PENAL, WEIGHT_ABSOLUTE_PENAL); - double penaltyQ = computeSlackPenalty(x, slackQStartIndex, numQEquations, WEIGHT_Q_PENAL, WEIGHT_ABSOLUTE_PENAL); - double penaltyV = computeSlackPenalty(x, slackVStartIndex, numVEquations, WEIGHT_V_PENAL, WEIGHT_ABSOLUTE_PENAL); - double totalPenalty = penaltyP + penaltyQ + penaltyV; - - LOGGER.info("==== Slack penalty details ===="); - LOGGER.info("Penalty P = {}", penaltyP); - LOGGER.info("Penalty Q = {}", penaltyQ); - LOGGER.info("Penalty V = {}", penaltyV); - LOGGER.info("Total penalty = {}", totalPenalty); - } - - /** - * Logs information like bus name and slack value for most significant slack variables. - * - * @param type The slack variable type. - * @param startIndex The start index of slack variables associated to the given type. - * @param count The maximum number of slack variables associated to the given type. - * @param x The variable values as returned by solver. - */ - private void logSlackValues(String type, int startIndex, int count, List x) { - final double sbase = 100.0; // Base power in MVA - - LOGGER.debug("==== Slack diagnostics for {} (p.u. and physical units) ====", type); - - for (int i = 0; i < count; i++) { - double sm = x.get(startIndex + 2 * i); - double sp = x.get(startIndex + 2 * i + 1); - double epsilon = sp - sm; - - // Get significant slack values above threshold - boolean shouldSkip = Math.abs(epsilon) <= knitroParameters.getSlackThreshold(); - String name = null; - String interpretation = null; - - if (!shouldSkip) { - name = getSlackVariableBusName(i, type); - - switch (type) { - case "P" -> interpretation = String.format("ΔP = %.4f p.u. (%.1f MW)", epsilon, epsilon * sbase); - case "Q" -> interpretation = String.format("ΔQ = %.4f p.u. (%.1f MVAr)", epsilon, epsilon * sbase); - case "V" -> { - var bus = network.getBusById(name); - if (bus == null) { - LOGGER.warn("Bus {} not found while logging V slack.", name); - shouldSkip = true; - } else { - interpretation = String.format("ΔV = %.4f p.u. (%.1f kV)", epsilon, epsilon * bus.getNominalV()); - } - } - default -> interpretation = "Unknown slack type"; - } - } - - if (shouldSkip) { - continue; - } - - String msg = String.format("Slack %s[ %s ] → Sm = %.4f, Sp = %.4f → %s", type, name, sm, sp, interpretation); - LOGGER.debug(msg); - } - } - - /** - * Finds the bus associated to a slack variable. - * - * @param index The index of the slack variable - * @param type The slack variable type. - * @return The id of the bus associated to the slack variable. - */ - private String getSlackVariableBusName(Integer index, String type) { - Set> equationSet = switch (type) { - case "P" -> pEquationLocalIds.entrySet(); - case "Q" -> qEquationLocalIds.entrySet(); - case "V" -> vEquationLocalIds.entrySet(); - default -> throw new IllegalStateException("Unexpected variable type: " + type); - }; - - Optional varIndexOptional = equationSet.stream() - .filter(entry -> index.equals(entry.getValue())) - .map(Map.Entry::getKey) - .findAny(); - - int varIndex; - if (varIndexOptional.isPresent()) { - varIndex = varIndexOptional.get(); - } else { - throw new PowsyblException("Variable index associated with slack variable " + type + " was not found"); - } - - LfBus bus = network.getBus(equationSystem.getIndex().getSortedEquationsToSolve().get(varIndex).getElementNum()); - - return bus.getId(); - } - - /** - * Calculates the total loss associated to a slack variable type - * - * @param x The variable values as returned by solver. - * @param startIndex The start index of slack variables associated to the given type. - * @param count The maximum number of slack variables associated to the given type. - * @param weight The weight inf front of the given slack variables terms - * @param lambda The coefficient of the linear terms in the objective function. - * @return The total penalty associated to the slack variables type. - */ - private double computeSlackPenalty(List x, int startIndex, int count, double weight, double lambda) { - double penalty = 0.0; - for (int i = 0; i < count; i++) { - double sm = x.get(startIndex + 2 * i); - double sp = x.get(startIndex + 2 * i + 1); - double diff = sp - sm; - penalty += weight * (diff * diff); // Quadratic terms - penalty += weight * lambda * (sp + sm); // Linear terms - } - return penalty; - } + public final class RelaxedKnitroProblem extends AbstractRelaxedKnitroProblem { - /** - * Optimization problem solving the open load flow equation system by minimizing constraint violations through relaxation. - */ - private final class RelaxedKnitroProblem extends AbstractKnitroProblem { + private RelaxedKnitroProblem(LfNetwork network, EquationSystem equationSystem, + TargetVector targetVector, JacobianMatrix jacobianMatrix, + KnitroSolverParameters knitroParameters, VoltageInitializer voltageInitializer) throws KNException { - /** - * Relaxed Knitro problem definition including: - * - initialization of variables (types, bounds, initial state) - * - definition of linear constraints - * - definition of non-linear constraints, evaluated in extended the callback function - * - definition of the extended Jacobian matrix passed to Knitro to solve the problem - * - definition of the objective function to be minimized (equation system violations) - */ - private RelaxedKnitroProblem( - LfNetwork network, - EquationSystem equationSystem, - TargetVector targetVector, - JacobianMatrix jacobianMatrix, - VoltageInitializer voltageInitializer, - KnitroSolverParameters parameters) throws KNException { + super(network, equationSystem, targetVector, jacobianMatrix, knitroParameters, numSlackVariables, 0); - super(network, equationSystem, targetVector, jacobianMatrix, parameters, numberOfVariables); - LOGGER.info("Defining {} variables", numberOfVariables); + LOGGER.info("Defining {} variables", numTotalVariables); // Initialize variables (base class handles LF variables, we customize for slack) - initializeVariables(voltageInitializer, numberOfVariables); + initializeVariables(voltageInitializer); LOGGER.info("Voltage initialization strategy: {}", voltageInitializer); // Set up the constraints of the relaxed optimization problem @@ -292,195 +77,7 @@ private RelaxedKnitroProblem( setJacobianMatrix(activeConstraints, nonlinearConstraintIndexes); // set the objective function of the optimization problem - // initialise lists to track quadratic objective function terms of the form: a * x1 * x2 - List quadRows = new ArrayList<>(); // list of indexes of the first variable x1 - List quadCols = new ArrayList<>(); // list of indexes of the second variable x2 - List quadCoefs = new ArrayList<>(); // list of indexes of the coefficient a - - // initialise lists to track linear objective function terms of the form: a * x - List linIndexes = new ArrayList<>(); // list of indexes of the variable x - List linCoefs = new ArrayList<>(); // list of indexes of the coefficient a - - // add slack penalty terms, for each slack type, of the form: (Sp - Sm)^2 = Sp^2 + Sm^2 - 2*Sp*Sm + linear terms from the absolute value - addSlackObjectiveTerms(numPEquations, slackPStartIndex, WEIGHT_P_PENAL, WEIGHT_ABSOLUTE_PENAL, quadRows, quadCols, quadCoefs, linIndexes, linCoefs); - addSlackObjectiveTerms(numQEquations, slackQStartIndex, WEIGHT_Q_PENAL, WEIGHT_ABSOLUTE_PENAL, quadRows, quadCols, quadCoefs, linIndexes, linCoefs); - addSlackObjectiveTerms(numVEquations, slackVStartIndex, WEIGHT_V_PENAL, WEIGHT_ABSOLUTE_PENAL, quadRows, quadCols, quadCoefs, linIndexes, linCoefs); - - setObjectiveQuadraticPart(quadRows, quadCols, quadCoefs); - setObjectiveLinearPart(linIndexes, linCoefs); - } - - /** - * Adds quadratic and linear terms related to slack variables to the objective function. - */ - private void addSlackObjectiveTerms(int numEquations, int slackStartIdx, double weight, double lambda, - List quadRows, List quadCols, List quadCoefs, - List linIndexes, List linCoefs) { - for (int i = 0; i < numEquations; i++) { - int idxSm = slackStartIdx + 2 * i; // negative slack variable index - int idxSp = slackStartIdx + 2 * i + 1; // positive slack variable index - - // Add quadratic terms: weight * (sp^2 + sm^2 - 2 * sp * sm) - - // add first quadratic term : weight * sp^2 - quadRows.add(idxSp); - quadCols.add(idxSp); - quadCoefs.add(weight); - - // add second quadratic term : weight * sm^2 - quadRows.add(idxSm); - quadCols.add(idxSm); - quadCoefs.add(weight); - - // add third quadratic term : weight * (- 2 * sp * sm) - quadRows.add(idxSp); - quadCols.add(idxSm); - quadCoefs.add(-2 * weight); - - // Add linear terms: weight * lambda * (sp + sm) - - // add first linear term : weight * lambda * sp - linIndexes.add(idxSp); - linCoefs.add(lambda * weight); - - // add second linear term : weight * lambda * sm - linIndexes.add(idxSm); - linCoefs.add(lambda * weight); - } - } - - @Override - protected void initializeCustomizedVariables(List lowerBounds, List upperBounds, - List initialValues, int numTotalVariables) { - // set a lower bound to slack variables (>= 0) - // initial values have already been set to 0 - for (int i = numberOfPowerFlowVariables; i < numTotalVariables; i++) { - lowerBounds.set(i, 0.0); - } - } - - @Override - protected void addAdditionalConstraintVariables(int equationId, AcEquationType equationType, - List varIndices, List coefficients) { - // Add slack variables if applicable - int slackBase = getSlackIndexBase(equationType, equationId); - if (slackBase >= 0) { - varIndices.add(slackBase); // Sm - varIndices.add(slackBase + 1); // Sp - coefficients.add(1.0); - coefficients.add(-1.0); - } - } - - @Override - protected void addAdditionalJacobianVariables(int constraintIndex, - Equation equation, - List variableIndices) { - AcEquationType equationType = equation.getType(); - // get slack variable local index (within its equation type) - int slackStart = switch (equationType) { - case BUS_TARGET_P -> pEquationLocalIds.getOrDefault(constraintIndex, -1); - case BUS_TARGET_Q -> qEquationLocalIds.getOrDefault(constraintIndex, -1); - case BUS_TARGET_V -> vEquationLocalIds.getOrDefault(constraintIndex, -1); - default -> -1; - }; - - if (slackStart >= 0) { - // get slack variable type starting index (within total variables' indexes) - int slackBaseIndex = switch (equationType) { - case BUS_TARGET_P -> slackPStartIndex; - case BUS_TARGET_Q -> slackQStartIndex; - case BUS_TARGET_V -> slackVStartIndex; - default -> throw new IllegalStateException("Unexpected constraint type: " + equationType); - }; - // get slack variables Sm and Sp indexes - variableIndices.add(slackBaseIndex + 2 * slackStart); // Sm - variableIndices.add(slackBaseIndex + 2 * slackStart + 1); // Sp - } - } - - /** - * Returns the base index of the slack variable associated with a given equation type and ID. - * - * @param equationType Type of the equation (P, Q, or V). - * @param equationId Index of the equation. - * @return Base index of the corresponding slack variable, or -1 if not applicable. - */ - private int getSlackIndexBase(AcEquationType equationType, int equationId) { - return switch (equationType) { - case BUS_TARGET_P -> pEquationLocalIds.getOrDefault(equationId, -1) >= 0 - ? slackPStartIndex + 2 * pEquationLocalIds.get(equationId) : -1; - case BUS_TARGET_Q -> qEquationLocalIds.getOrDefault(equationId, -1) >= 0 - ? slackQStartIndex + 2 * qEquationLocalIds.get(equationId) : -1; - case BUS_TARGET_V -> vEquationLocalIds.getOrDefault(equationId, -1) >= 0 - ? slackVStartIndex + 2 * vEquationLocalIds.get(equationId) : -1; - default -> -1; - }; - } - - @Override - protected KNEvalGACallback createGradientCallback(JacobianMatrix jacobianMatrix, - List listNonZerosCtsDense, List listNonZerosVarsDense, - List listNonZerosCtsSparse, List listNonZerosVarsSparse) { - - return new RelaxedCallbackEvalG(jacobianMatrix, listNonZerosCtsDense, listNonZerosVarsDense, - listNonZerosCtsSparse, listNonZerosVarsSparse, network, - equationSystem, knitroParameters, numberOfPowerFlowVariables); - } - - /** - * Callback used by Knitro to evaluate the non-linear parts of the objective and constraint functions. - */ - private static final class RelaxedCallbackEvalFC extends KnitroCallbacks.BaseCallbackEvalFC { - - private final RelaxedKnitroProblem problemInstance; - - private RelaxedCallbackEvalFC(RelaxedKnitroProblem problemInstance, - List> sortedEquationsToSolve, - List nonLinearConstraintIds) { - super(sortedEquationsToSolve, nonLinearConstraintIds); - this.problemInstance = problemInstance; - } - - @Override - protected double addModificationOfNonLinearConstraints(int equationId, AcEquationType equationType, - List x) { - int slackIndexBase = problemInstance.getSlackIndexBase(equationType, equationId); - if (slackIndexBase >= 0) { - double sm = x.get(slackIndexBase); // negative slack - double sp = x.get(slackIndexBase + 1); // positive slack - return sp - sm; // add slack contribution - } - return 0; - } - } - - /** - * Callback used by Knitro to evaluate the gradient (Jacobian matrix) of the constraints. - * Only constraints (no objective) are handled here. - */ - private static final class RelaxedCallbackEvalG extends KnitroCallbacks.BaseCallbackEvalG { - - private RelaxedCallbackEvalG(JacobianMatrix jacobianMatrix, - List denseConstraintIndices, List denseVariableIndices, - List sparseConstraintIndices, List sparseVariableIndices, - LfNetwork network, EquationSystem equationSystem, - KnitroSolverParameters knitroParameters, int numLFVariables) { - - super(jacobianMatrix, denseConstraintIndices, denseVariableIndices, sparseConstraintIndices, sparseVariableIndices, - network, equationSystem, knitroParameters, numLFVariables); - } - - @Override - protected double computeModifiedJacobianValue(int variableIndex, int constraintIndex) { - if (((variableIndex - numLFVariables) & 1) == 0) { - // set Jacobian entry to -1.0 if slack variable is Sm - return -1.0; - } else { - // set Jacobian entry to +1.0 if slack variable is Sp - return 1.0; - } - } + addObjectiveFunction(numPEquations, slackPStartIndex, numQEquations, slackQStartIndex, numVEquations, slackVStartIndex); } } } diff --git a/src/main/java/com/powsybl/openloadflow/knitro/solver/UseReactiveLimitsKnitroSolver.java b/src/main/java/com/powsybl/openloadflow/knitro/solver/UseReactiveLimitsKnitroSolver.java new file mode 100644 index 00000000..5bc180da --- /dev/null +++ b/src/main/java/com/powsybl/openloadflow/knitro/solver/UseReactiveLimitsKnitroSolver.java @@ -0,0 +1,637 @@ +/** + * Copyright (c) 2025, Artelys (http://www.artelys.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/. + * SPDX-License-Identifier: MPL-2.0 + */ +package com.powsybl.openloadflow.knitro.solver; + +import com.artelys.knitro.api.*; +import com.artelys.knitro.api.callbacks.KNEvalGACallback; +import com.powsybl.commons.PowsyblException; +import com.powsybl.openloadflow.ac.equations.AcEquationType; +import com.powsybl.openloadflow.ac.equations.AcVariableType; +import com.powsybl.openloadflow.equations.*; +import com.powsybl.openloadflow.network.*; +import com.powsybl.openloadflow.network.util.VoltageInitializer; +import org.slf4j.Logger; +import org.slf4j.LoggerFactory; + +import java.util.*; +import java.util.stream.Stream; + +import static com.powsybl.openloadflow.ac.equations.AcEquationType.*; + +/** + * Relaxed Knitro solver, solving the open load flow equation system by minimizing constraint violations through relaxation, + * and taking into account generator reactive limits. + * The consideration of limits is ensured by adding complementarity constraints in the problem formulation. + * The solver must therefore extend the open load flow equation system. + * + * It should be noted that complementarity constraints fix the reactive limit on one side if it is violated, + * and relax the voltage setpoint through the addition of auxiliary variables. + * It should be noted that this relaxation is not penalized, as it corresponds to a PV/PQ transition of the bus. + * Note that voltage setpoints remain relaxed through the addition of a slack variable. + * + * @author Pierre Arvy {@literal } + * @author Yoann Anezin {@literal } + */ +public class UseReactiveLimitsKnitroSolver extends AbstractRelaxedKnitroSolver { + + private static final Logger LOGGER = LoggerFactory.getLogger(UseReactiveLimitsKnitroSolver.class); + + // threshold to identify well-defined generator reactive power limits + private static final double MAX_REASONABLE_REACTIVE_LIMIT = 1e15; + + // number of buses for which the reactive limits are well-defined + // for each of these equations, some complementarity constraints / variables will be added to Knitro optimization problem, + // to model the PV/PQ switches of generators + private final int numBusesWithFiniteQLimits; + + // Starting indices for slack variables in the variable vector + private final int compVarStartIndex; + + // Mappings from global equation indices to local indices by equation type + private final Map elemNumControlledControllerBus; + private static final Map> INDEQUNACTIVEQ = new LinkedHashMap<>(); + + // reactive limit equations to add to the constraints to model the PV/PQ switches of generators + private final List> equationsQBusV; + private final List busesNumWithReactiveLimitEquationsToAdd; + + // mappings from global equation indices to local indices for the voltage target equations to duplicate + private final Map vSuppEquationLocalIds; + + public UseReactiveLimitsKnitroSolver(LfNetwork network, KnitroSolverParameters knitroParameters, EquationSystem equationSystem, + JacobianMatrix j, TargetVector targetVector, + EquationVector equationVector, boolean detailedReport) { + super(network, knitroParameters, equationSystem, j, targetVector, equationVector, detailedReport); + + // the optimization problem modeled here duplicates and modifies V equations of open load flow equations system, to model + // the voltage change due to a PV/PQ switch of a generator + // to do this, it also adds two Q equations on each PV bus, fixing the setpoint to a limit, if the reactive power exceeds this limit + this.compVarStartIndex = slackVStartIndex + 2 * numVEquations; + + // to define a coherent system, first we need to collect the Q equations for which reactive bounds are well-defined + this.elemNumControlledControllerBus = new LinkedHashMap<>(); + + // first, we retrieve the V equations from the open load flow system, because reactive power limits are added + // for generators that control voltage + List> olfSortedEquations = equationSystem.getIndex().getSortedEquationsToSolve(); + List busesWithVoltageTargetEquation = olfSortedEquations.stream() + .filter(e -> e.getType() == AcEquationType.BUS_TARGET_V) + .map(e -> e.getTerms().getFirst().getElementNum()) + .toList(); + + // contains the reactive equations that must be added to the optimization problem constraints + List> reactiveEquationsToAdd = new ArrayList<>(); + // contains the buses for which reactive limit equation must be added to the constraints + this.busesNumWithReactiveLimitEquationsToAdd = new ArrayList<>(); + for (int elementNum : busesWithVoltageTargetEquation) { + LfBus controlledBus = network.getBuses().get(elementNum); + + // the reactive limits is only supported for buses whose voltage is controlled by a generator + controlledBus.getGeneratorVoltageControl().ifPresent( + generatorVoltageControl -> { + LfBus controllerBus = generatorVoltageControl.getControllerElements().getFirst(); + + // only buses with well-defined reactive power limits are considered + // otherwise, we cannot fix the reactive power via complementarity, in the optimization problem modeled here + if (Double.isFinite(controllerBus.getMaxQ()) && Math.abs(controllerBus.getMaxQ()) < MAX_REASONABLE_REACTIVE_LIMIT + && Double.isFinite(controllerBus.getMinQ()) && Math.abs(controllerBus.getMinQ()) < MAX_REASONABLE_REACTIVE_LIMIT) { + + // retrieve the reactive power balance equation of the bus that controls voltage + // the equation is inactive in the open load flow equation system, but it will be used to evaluate the reactive power balance of the generator + // in the optimization problem constraints + List> controllerBusEquations = equationSystem.getEquations(ElementType.BUS, controllerBus.getNum()); + Equation reactiveEquationToAdd = controllerBusEquations.stream() + .filter(e -> e.getType() == BUS_TARGET_Q) + .toList() + .getFirst(); + + // add the equations to model in the corresponding lists + reactiveEquationsToAdd.add(reactiveEquationToAdd); + elemNumControlledControllerBus.put(controllerBus.getNum(), controlledBus.getNum()); // link between controller and controlled bus + busesNumWithReactiveLimitEquationsToAdd.add(controllerBus.getNum()); + } + } + ); + } + + // this will be used to iterate on the equations to add in the system + this.numBusesWithFiniteQLimits = reactiveEquationsToAdd.size(); + + // duplication to later model the equations defining the bLow and bUp variables, equal + // to the equations modeling the fixing of reactive power limits + this.equationsQBusV = Stream.concat(reactiveEquationsToAdd.stream(), reactiveEquationsToAdd.stream()).toList(); + + // map added voltage target equations to local indices + this.vSuppEquationLocalIds = new HashMap<>(); + int[] vSuppCounter = {0}; + for (int i = 0; i < olfSortedEquations.size(); i++) { + final int equationIndex = i; + Equation equation = olfSortedEquations.get(i); + AcEquationType type = equation.getType(); + if (Objects.requireNonNull(type) == BUS_TARGET_V) { + // add a vSup equation + equation.getElement(network).ifPresent( + element -> { + LfBus controlledBus = network.getBuses().get(element.getNum()); + // supports only voltage control of generators + controlledBus.getGeneratorVoltageControl().ifPresent( + generatorVoltageControl -> { + LfBus controllerBus = generatorVoltageControl.getControllerElements().getFirst(); + if (busesNumWithReactiveLimitEquationsToAdd.contains(controllerBus.getNum())) { + vSuppEquationLocalIds.put(equationIndex, vSuppCounter[0]++); + } + } + ); + } + ); + } + } + } + + /** + * Returns the name of the solver. + */ + @Override + public String getName() { + return "Knitro Generator Reactive Limits Solver"; + } + + @Override + protected KNProblem createKnitroProblem(VoltageInitializer voltageInitializer) { + try { + return new UseReactiveLimitsKnitroProblem(network, equationSystem, targetVector, j, knitroParameters, voltageInitializer); + } catch (KNException e) { + throw new PowsyblException("Failed to create Knitro problem modeling generator reactive limits", e); + } + } + + @Override + protected void processSolution(KNSolver solver, KNSolution solution, KNProblem problemInstance) { + // log solution, slacks and penalties + super.processSolution(solver, solution, problemInstance); + + // log the switches computed by the optimization + logSwitches(solution, compVarStartIndex, numBusesWithFiniteQLimits); + } + + /** + * Inform all switches PV -> PQ done in the solution found + * + * @param solution Solution returned by the optimization + * @param startIndex first index of complementarity constraints variables in x + * @param count number of b_low / b_up different variables + */ + private void logSwitches(KNSolution solution, int startIndex, int count) { + LOGGER.info("=== Switches Done==="); + List x = solution.getX(); // solution returned by the optimization + double eps = knitroParameters.getSlackThreshold(); + + for (int i = 0; i < count; i++) { + double vInf = x.get(startIndex + 5 * i); + double vSup = x.get(startIndex + 5 * i + 1); + double bLow = x.get(startIndex + 5 * i + 2); + double bUp = x.get(startIndex + 5 * i + 3); + + int controllerBusNum = busesNumWithReactiveLimitEquationsToAdd.get(i); + // get the bus ID from the network + String busId = network.getBuses().get(controllerBusNum).getId(); + + // switch to lower bound case: bLow is null and the auxiliary variable is not + if (Math.abs(bLow) < eps && (vInf > eps || vSup > eps)) { + LOGGER.info("Switch PV -> PQ on bus {}, Q set at Qmin", busId); + + // switch to upper bound case: bUp is null and the auxiliary variable is not + } else if (Math.abs(bUp) < eps && (vInf > eps || vSup > eps)) { + LOGGER.info("Switch PV -> PQ on bus {}, Q set at Qmax", busId); + } + } + } + + public final class UseReactiveLimitsKnitroProblem extends AbstractRelaxedKnitroProblem { + + // objects to retrieve the complete form of the equation system that we seek to satisfy as constraints + // of the optimization problem + // to model generator reactive limits, additional equations and therefore right-hand sides are necessary, and these + // are not considered in the open load flow objects + private List> completeEquationsToSolve; + private List completeTargetVector; + + /** + * Generator reactive limits Knitro problem definition including: + * - Initialization of customized variables (complmentarity auxiliary variables, and those specific to the Knitro modeling) + * - Definition of linear and non-linear constraints + * - Jacobian matrix setup for Knitro + */ + private UseReactiveLimitsKnitroProblem(LfNetwork network, EquationSystem equationSystem, + TargetVector targetVector, JacobianMatrix jacobianMatrix, + KnitroSolverParameters parameters, VoltageInitializer voltageInitializer) throws KNException { + // initialize optimization problem + + // for each complementary equation to add, 3 new variables are used on V equations and 2 on Q equations + // these "auxiliary variables" allow modeling the relaxation of set points or limits + // the total number of variables includes these variables, as well as those from the open load flow system + // and from the system relaxation (slack variables) + + // the total number of equations equals the number of constraints in the open load flow system, + // to which 3 equations are added for each voltage control equation whose associated generator + // has finite reactive power limits + super(network, equationSystem, targetVector, jacobianMatrix, parameters, numSlackVariables + numBusesWithFiniteQLimits * 5, 3 * numBusesWithFiniteQLimits); + + LOGGER.info("Defining {} variables", numTotalVariables); + + // Set up the constraints of the optimization problem + setupConstraints(); + + // Initialize variables + initializeVariables(voltageInitializer); + LOGGER.info("Initialization of variables : type of initialization {}", voltageInitializer); + + // specify the callback to evaluate the jacobian + setObjEvalCallback(new UseReactiveLimitsCallbackEvalFC(this, completeEquationsToSolve, nonlinearConstraintIndexes)); + + // set the pattern of the jacobian + setJacobianMatrix(completeEquationsToSolve, nonlinearConstraintIndexes); + + // set the objective function of the optimization problem + addObjectiveFunction(numPEquations, slackPStartIndex, numQEquations, slackQStartIndex, numVEquations, slackVStartIndex); + } + + @Override + protected void initializeCustomizedVariables(List lowerBounds, List upperBounds, List initialValues) { + // initialize slack variables + super.initializeCustomizedVariables(lowerBounds, upperBounds, initialValues); + + // initialize auxiliary variables in complementary constraints + for (int i = 0; i < numBusesWithFiniteQLimits; i++) { + lowerBounds.set(compVarStartIndex + 5 * i, 0.0); + lowerBounds.set(compVarStartIndex + 5 * i + 1, 0.0); + + lowerBounds.set(compVarStartIndex + 5 * i + 2, 0.0); + lowerBounds.set(compVarStartIndex + 5 * i + 3, 0.0); + lowerBounds.set(compVarStartIndex + 5 * i + 4, 0.0); + + initialValues.set(compVarStartIndex + 5 * i + 2, 1.0); + initialValues.set(compVarStartIndex + 5 * i + 3, 1.0); + } + } + + @Override + protected void setUpScalingFactors() throws KNException { + List scalingFactors = new ArrayList<>(Collections.nCopies(numTotalVariables, 1.0)); + List scalingCenters = new ArrayList<>(Collections.nCopies(numTotalVariables, 0.0)); + for (int i = numberOfPowerFlowVariables; i < slackVStartIndex; i++) { + scalingFactors.set(i, 1e-2); + } + + ArrayList list = new ArrayList<>(numTotalVariables); + for (int i = 0; i < numTotalVariables; i++) { + list.add(i); + } + + setVarScaleFactors(new KNSparseVector<>(list, scalingFactors)); + setVarScaleCenters(new KNSparseVector<>(list, scalingCenters)); + } + + @Override + protected void setupConstraints() throws KNException { + activeConstraints = equationSystem.getIndex().getSortedEquationsToSolve(); + + // this problem models complementarity constraints to model the PV/PQ switching of buses + // to do this, constraints are added to the system compared to the OLF system, and therefore right-hand sides are also added + // it is therefore not possible to use the system and right-hand side directly, which is why we initialize the following objects + // these include all the elements of the OLF system / rhs, and all those that are added (for complementarity constraints) + completeEquationsToSolve = new ArrayList<>(activeConstraints); // contains all equations of the final system to be solved + completeTargetVector = new ArrayList<>(Arrays.stream(targetVector.getArray()).boxed().toList()); // contains all the target of the system to be solved + + // create solver utils here to only create one + NonLinearExternalSolverUtils solverUtils = new NonLinearExternalSolverUtils(); + + // add linear constraints and fill the list of non-linear constraints + addLinearConstraints(activeConstraints, solverUtils); + + // number of variables from the OLF system and whose voltage control equations have been duplicated + // for cases where reactive power limits are well-defined + int intermediateNumberOfActiveEquations = completeEquationsToSolve.size(); + + // add the constraints that allow considering reactive power limits + completeEquationsToSolve.addAll(equationsQBusV); + for (int equationId = intermediateNumberOfActiveEquations; equationId < completeEquationsToSolve.size(); equationId++) { + Equation equation = completeEquationsToSolve.get(equationId); + LfBus controllerBus = network.getBus(equation.getElementNum()); + if (equationId - intermediateNumberOfActiveEquations < equationsQBusV.size() / 2) { + completeTargetVector.add(controllerBus.getMinQ() - controllerBus.getLoadTargetQ()); // bLow target + } else { + completeTargetVector.add(controllerBus.getMaxQ() - controllerBus.getLoadTargetQ()); // bUp target + } + // these equations are non-linear since reactive power flows appear in the balance + nonlinearConstraintIndexes.add(equationId); + INDEQUNACTIVEQ.put(equationId, equation); + } + + // this includes the number of equations from the OLF system, + // but also all those added to model complementarity + int totalNumConstraints = completeEquationsToSolve.size(); + LOGGER.info("Defining {} active constraints", totalNumConstraints); + + // pass to Knitro the indexes of non-linear constraints, that will be evaluated in the callback function + // NOTE: in the non-linear constraints, there are also here specific elements for the constraints that we have added + setMainCallbackCstIndexes(nonlinearConstraintIndexes); + + // right hand side (targets), specific to the complete problem we are trying to solve + setConEqBnds(completeTargetVector); + + // declaration of complementarity constraints in Knitro + addComplementaryConstraints(); + } + + /** + * Adds a single constraint to the Knitro problem. + * Linear constraints are directly encoded; non-linear ones are delegated to the callback. + * + * @param equationId Index of the equation in the list. + */ + @Override + protected void addConstraint(int equationId, List> sortedEquationsToSolve, + NonLinearExternalSolverUtils solverUtils) { + + Equation equation = sortedEquationsToSolve.get(equationId); + AcEquationType equationType = equation.getType(); + List> terms = equation.getTerms(); + + if (equationType == BUS_TARGET_V) { + equation.getElement(network).ifPresent(lfElement -> { + LfBus controlledBus = network.getBuses().get(lfElement.getNum()); + controlledBus.getGeneratorVoltageControl().ifPresent(generatorVoltageControl -> { + LfBus controllerBus = generatorVoltageControl.getControllerElements().getFirst(); + + // boolean to indicate if the V equation will be duplicated + boolean addComplementarityConstraintsVariable = busesNumWithReactiveLimitEquationsToAdd.contains(controllerBus.getNum()); + + if (NonLinearExternalSolverUtils.isLinear(equationType, terms)) { + try { + // Extract linear constraint components + var linearConstraint = solverUtils.getLinearConstraint(equationType, terms); + List varVInfIndices = new ArrayList<>(linearConstraint.listIdVar()); + List coefficientsVInf = new ArrayList<>(linearConstraint.listCoef()); + + // To add complementarity conditions, Knitro requires that they be written as two variables + // that complement each other. That is why we are introducing new variables that will play this role. + // We call them V_inf, V_sup, b_low and b_up. The two lasts appear in non-linear constraints + // Equations on V are duplicated, we add V_inf to one and V_sup to the other. + // We are also adding a variable V_aux that allows us to perform the PV -> PQ switch. + + // vInf equation, from the OLF system equation + // this equation serves to relax the voltage constraint when the maximum reactive power limit is reached + if (addComplementarityConstraintsVariable) { + int compVarBaseIndex = getComplementarityVarBaseIndex(equationId); + varVInfIndices.add(compVarBaseIndex); // vInf + varVInfIndices.add(compVarBaseIndex + 4); // vAux + coefficientsVInf.add(1.0); + coefficientsVInf.add(-1.0); + } + + // add slack variables if applicable + // NOTE: slacks can be added to relax the equation, even if there is no complementarity for this equation + addAdditionalConstraintVariables(equationId, equationType, varVInfIndices, coefficientsVInf); + + for (int i = 0; i < varVInfIndices.size(); i++) { + this.addConstraintLinearPart(equationId, varVInfIndices.get(i), coefficientsVInf.get(i)); + } + + LOGGER.trace("Added linear constraint #{} of type {}", equationId, equationType); + + // vSup equation, added only if it is a complementarity constraint + // this equation serves to relax the voltage constraint when the minimum reactive power limit is reached + if (addComplementarityConstraintsVariable) { + List varVSupIndices = new ArrayList<>(linearConstraint.listIdVar()); + List coefficientsVSup = new ArrayList<>(linearConstraint.listCoef()); + + // add complementarity constraints' variables + int compVarBaseIndex = getComplementarityVarBaseIndex(equationId); + varVSupIndices.add(compVarBaseIndex + 1); // V_sup + varVSupIndices.add(compVarBaseIndex + 4); // V_aux + coefficientsVSup.add(-1.0); + coefficientsVSup.add(1.0); + + // add slack variables if applicable + addAdditionalConstraintVariables(equationId, equationType, varVSupIndices, coefficientsVSup); + + for (int i = 0; i < varVSupIndices.size(); i++) { + int equationIndex = sortedEquationsToSolve.size() + vSuppEquationLocalIds.get(equationId); + this.addConstraintLinearPart(equationIndex, varVSupIndices.get(i), coefficientsVSup.get(i)); + } + } + } catch (UnsupportedOperationException e) { + throw new PowsyblException("Failed to process linear constraint for equation #" + equationId, e); + } + + // if the V equation is not linear, we add the index to non-linear equations (as well as the duplicate vSup) + } else { + nonlinearConstraintIndexes.add(equationId); + nonlinearConstraintIndexes.add(sortedEquationsToSolve.size() + vSuppEquationLocalIds.get(equationId)); + } + + // the duplicated voltage equations have been specified to the Knitro problem in the previous lines + // in order to maintain consistency with the modeled equation system, as well as the rhs, we must also update the corresponding objects + // in particular, we must add the equation for vSup to the list of modeled equations and to the rhs + if (addComplementarityConstraintsVariable) { + completeEquationsToSolve.add(equation); + completeTargetVector.add(Arrays.stream(targetVector.getArray()).boxed().toList().get(equationId)); // we apply the same voltage + } + }); + }); + // for other type of equations, the constraint can be added as usual + } else { + super.addConstraint(equationId, sortedEquationsToSolve, solverUtils); + } + } + + /** + * Adds the complementarity constraints to the Knitro problem definition. + */ + private void addComplementaryConstraints() throws KNException { + List listTypeVar = new ArrayList<>(Collections.nCopies(2 * numBusesWithFiniteQLimits, KNConstants.KN_CCTYPE_VARVAR)); + List bVarList = new ArrayList<>(); // bUp and bLow + List vInfSuppList = new ArrayList<>(); // vInf and vSup + + for (int i = 0; i < numBusesWithFiniteQLimits; i++) { + vInfSuppList.add(compVarStartIndex + 5 * i); // vInf + vInfSuppList.add(compVarStartIndex + 5 * i + 1); // vSup + bVarList.add(compVarStartIndex + 5 * i + 3); // bUp + bVarList.add(compVarStartIndex + 5 * i + 2); // bLow + } + + // add the complementary conditions + // 0 <= bLow perp vSup >= 0 and 0 <= bUp perp vInf >= 0 + setCompConstraintsTypes(listTypeVar); + setCompConstraintsParts(bVarList, vInfSuppList); + } + + private int getComplementarityVarBaseIndex(int equationId) { + return compVarStartIndex + 5 * vSuppEquationLocalIds.get(equationId); + } + + private int getElemNumControlledBus(int elemNum) { + return elemNumControlledControllerBus.get(elemNum); + } + + @Override + protected KNEvalGACallback createGradientCallback(JacobianMatrix jacobianMatrix, + List listNonZerosCtsDense, List listNonZerosVarsDense, + List listNonZerosCtsSparse, List listNonZerosVarsSparse) { + + return new UseReactiveLimitsCallbackEvalG(jacobianMatrix, listNonZerosCtsDense, listNonZerosVarsDense, + listNonZerosCtsSparse, listNonZerosVarsSparse, network, equationSystem, + knitroParameters, numberOfPowerFlowVariables); + } + + @Override + protected void addAdditionalJacobianVariables(int constraintIndex, + Equation equation, + List variableIndices) { + super.addAdditionalJacobianVariables(constraintIndex, equation, variableIndices); + int numberLFEq = completeEquationsToSolve.size() - 3 * vSuppEquationLocalIds.size(); + AcEquationType equationType = equation.getType(); + // Add complementarity constraints' variables if the constraint type has them + int compVarStart; + // Case of inactive Q equations, we take the V equation associated to it to get the right index used to order the equation system + if (equationType == BUS_TARGET_Q && !equation.isActive()) { + int elemNumControlledBus = elemNumControlledControllerBus.get(equation.getElementNum()); // Controller bus + List> listEqControlledBus = equationSystem // Equations of the Controller bus + .getEquations(ElementType.BUS, elemNumControlledBus); + Equation eqVControlledBus = listEqControlledBus.stream() // Take the one on V + .filter(e -> e.getType() == BUS_TARGET_V).toList().getFirst(); + int indexEqVAssociated = equationSystem.getIndex().getSortedEquationsToSolve().indexOf(eqVControlledBus); // Find the index of the V equation associated + + compVarStart = vSuppEquationLocalIds.get(indexEqVAssociated); + if (constraintIndex < numberLFEq + 2 * vSuppEquationLocalIds.size()) { + variableIndices.add(compVarStartIndex + 5 * compVarStart + 2); // b_low + } else { + variableIndices.add(compVarStartIndex + 5 * compVarStart + 3); // b_up + } + } + } + + /** + * Callback used by Knitro to evaluate the non-linear parts of the objective and constraint functions. + */ + private static final class UseReactiveLimitsCallbackEvalFC extends RelaxedCallbackEvalFC { + + private final UseReactiveLimitsKnitroProblem problemInstance; + + private UseReactiveLimitsCallbackEvalFC(UseReactiveLimitsKnitroProblem problemInstance, + List> sortedEquationsToSolve, + List nonLinearConstraintIds) { + super(problemInstance, sortedEquationsToSolve, nonLinearConstraintIds); + this.problemInstance = problemInstance; + } + + @Override + protected double addModificationOfNonLinearConstraints(int equationId, AcEquationType equationType, + List x) { + double constraintValue = 0; + Equation equation = sortedEquationsToSolve.get(equationId); + // if the equation is active, then it is treated as an open load flow constraint + if (equation.isActive()) { + // add relaxation if necessary + constraintValue += super.addModificationOfNonLinearConstraints(equationId, equationType, x); + + // otherwise, these are bLow and bUp constraints and the term must be added + } else { + int numNotActiveEquations = sortedEquationsToSolve.stream().filter( + e -> !e.isActive()).toList().size(); + int elemNum = equation.getElementNum(); + int elemNumControlledBus = problemInstance.getElemNumControlledBus(elemNum); + List> controlledBusEquations = sortedEquationsToSolve.stream() + .filter(e -> e.getElementNum() == elemNumControlledBus).toList(); + // the voltage set point equation + Equation equationV = controlledBusEquations.stream().filter(e -> e.getType() == BUS_TARGET_V).toList().getFirst(); + int equationVId = sortedEquationsToSolve.indexOf(equationV); //Index of V equation + int compVarBaseIndex = problemInstance.getComplementarityVarBaseIndex(equationVId); + // bLow constraint + if (equationId - sortedEquationsToSolve.size() + numNotActiveEquations / 2 < 0) { + double bLow = x.get(compVarBaseIndex + 2); + constraintValue -= bLow; + // bUp Constraint + } else { + double bUp = x.get(compVarBaseIndex + 3); + constraintValue += bUp; + } + } + + return constraintValue; + } + } + + /** + * Callback used by Knitro to evaluate the gradient (Jacobian matrix) of the constraints. + * Only constraints (no objective) are handled here. + */ + private static final class UseReactiveLimitsCallbackEvalG extends RelaxedCallbackEvalG { + + private final int numLFVar; + private final int numVEq; + private final int numPQEq; + + private final Map> indRowVariable = new LinkedHashMap<>(); + + private UseReactiveLimitsCallbackEvalG(JacobianMatrix jacobianMatrix, + List denseConstraintIndices, List denseVariableIndices, List sparseConstraintIndices, + List sparseVariableIndices, LfNetwork network, EquationSystem equationSystem, + KnitroSolverParameters knitroParameters, int numLFVariables) { + + super(jacobianMatrix, denseConstraintIndices, denseVariableIndices, sparseConstraintIndices, sparseVariableIndices, + network, equationSystem, knitroParameters, numLFVariables); + + this.numLFVar = equationSystem.getIndex().getSortedVariablesToFind().size(); + + this.numVEq = equationSystem.getIndex().getSortedEquationsToSolve().stream().filter( + e -> e.getType() == BUS_TARGET_V).toList().size(); + + this.numPQEq = equationSystem.getIndex().getSortedEquationsToSolve().stream().filter( + e -> e.getType() == BUS_TARGET_Q || e.getType() == BUS_TARGET_P).toList().size(); + + List> listVar = equationSystem.getVariableSet().getVariables().stream().toList(); + for (Variable variable : listVar) { + indRowVariable.put(variable.getRow(), variable); + } + } + + @Override + protected double computeAddedConstraintJacobianValue(int variableIndex, int constraintIndex) { + // Case of Unactivated Q equations + // If var is a LF variable : derivate non-activated equations + double value = 0.0; + Equation equation = INDEQUNACTIVEQ.get(constraintIndex); + if (variableIndex < numLFVar) { + // add non-linear terms + for (Map.Entry, List>> e : equation.getTermsByVariable().entrySet()) { + for (EquationTerm term : e.getValue()) { + Variable v = e.getKey(); + if (indRowVariable.get(variableIndex) == v) { + value += term.isActive() ? term.der(v) : 0; + } + } + } + } + // Check if var is a b_low or b_up var + if (variableIndex >= numLFVar + 2 * (numPQEq + numVEq)) { + int rest = (variableIndex - numLFVar - 2 * (numPQEq + numVEq)) % 5; + if (rest == 2) { + // set Jacobian entry to -1.0 if variable is b_low + value = -1.0; + } else if (rest == 3) { + // set Jacobian entry to 1.0 if variable is b_up + value = 1.0; + } + } + return value; + } + } + } +} diff --git a/src/test/java/com/powsybl/openloadflow/knitro/solver/UseReactiveLimitsKnitroSolverTest.java b/src/test/java/com/powsybl/openloadflow/knitro/solver/UseReactiveLimitsKnitroSolverTest.java new file mode 100644 index 00000000..12c7d169 --- /dev/null +++ b/src/test/java/com/powsybl/openloadflow/knitro/solver/UseReactiveLimitsKnitroSolverTest.java @@ -0,0 +1,307 @@ +/** + * Copyright (c) 2025, Artelys (http://www.artelys.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/. + * SPDX-License-Identifier: MPL-2.0 + */ +package com.powsybl.openloadflow.knitro.solver; + +import com.powsybl.iidm.network.Network; +import com.powsybl.ieeecdf.converter.IeeeCdfNetworkFactory; +import com.powsybl.iidm.network.*; +import com.powsybl.iidm.network.test.EurostagTutorialExample1Factory; +import com.powsybl.loadflow.LoadFlow; +import com.powsybl.loadflow.LoadFlowParameters; +import com.powsybl.loadflow.LoadFlowResult; +import com.powsybl.math.matrix.DenseMatrixFactory; +import com.powsybl.openloadflow.OpenLoadFlowParameters; +import com.powsybl.openloadflow.OpenLoadFlowProvider; +import com.powsybl.openloadflow.network.EurostagFactory; +import org.junit.jupiter.api.BeforeEach; +import org.junit.jupiter.api.Test; + +import java.util.concurrent.CompletionException; + +import static com.powsybl.openloadflow.util.LoadFlowAssert.assertReactivePowerEquals; +import static org.junit.jupiter.api.Assertions.*; + +/** + * @author Pierre Arvy {@literal } + * @author Yoann Anezin {@literal } + * @author Amine Makhen {@literal } + */ +class UseReactiveLimitsKnitroSolverTest { + private LoadFlow.Runner loadFlowRunner; + private LoadFlowParameters parameters; + private KnitroLoadFlowParameters knitroParams; + + @BeforeEach + void setUp() { + loadFlowRunner = new LoadFlow.Runner(new OpenLoadFlowProvider(new DenseMatrixFactory())); + + parameters = new LoadFlowParameters() + // no need to activate reactive limits outer loop as this is directly supported by use reactive limits knitro solver + .setUseReactiveLimits(false) + .setDistributedSlack(false); + + OpenLoadFlowParameters.create(parameters) + .setAcSolverType(KnitroSolverFactory.NAME) + .setVoltageInitModeOverride(OpenLoadFlowParameters.VoltageInitModeOverride.FULL_VOLTAGE); + + knitroParams = new KnitroLoadFlowParameters() + .setKnitroSolverType(KnitroSolverParameters.SolverType.USE_REACTIVE_LIMITS) + .setGradientComputationMode(2) + .setThreadNumber(1); + + parameters.addExtension(KnitroLoadFlowParameters.class, knitroParams); + } + + private static Network modifiedEurostagFactory(double qLow, double qUp) { + Network network = EurostagFactory.fix(EurostagTutorialExample1Factory.create()); + + // access to already created equipments + Generator gen = network.getGenerator("GEN"); + Substation p1 = network.getSubstation("P1"); + + // reduce GEN reactive range + gen.newMinMaxReactiveLimits() + .setMinQ(0) + .setMaxQ(280) + .add(); + + // create a new generator GEN2 + VoltageLevel vlgen2 = p1.newVoltageLevel() + .setId("VLGEN2") + .setNominalV(24.0) + .setTopologyKind(TopologyKind.BUS_BREAKER) + .add(); + + vlgen2.getBusBreakerView() + .newBus() + .setId("NGEN2") + .add(); + + vlgen2.newGenerator() + .setId("GEN2") + .setBus("NGEN2") + .setConnectableBus("NGEN2") + .setMinP(-9999.99) + .setMaxP(9999.99) + .setVoltageRegulatorOn(true) + .setTargetV(24.5) + .setTargetP(100) + .add(); + + // update added generator reactive limits + network.getGenerator("GEN2") + .newMinMaxReactiveLimits() + .setMinQ(qLow) + .setMaxQ(qUp) + .add(); + + p1.newTwoWindingsTransformer() + .setId("NGEN2_NHV1") + .setBus1("NGEN2") + .setConnectableBus1("NGEN2") + .setRatedU1(24.0) + .setBus2("NHV1") + .setConnectableBus2("NHV1") + .setRatedU2(400.0) + .setR(0.19253333333333333) + .setX(8.019911489428772) + .add(); + + // fix active power balance + network.getLoad("LOAD") + .setP0(699.838); + + return network; + } + + /** + * Modification of the network to urge a PV-PQ switch and set the new PQ bus to the lower bound on reactive power + */ + @Test + void testReacLimEurostagQlow() { + Network network = modifiedEurostagFactory(250, 300); + + // verify convergence using finite differences + LoadFlowResult result = loadFlowRunner.run(network, parameters); + assertTrue(result.isFullyConverged()); + + Generator gen2 = network.getGenerator("GEN2"); + TwoWindingsTransformer ngen2Nhv1 = network.getTwoWindingsTransformer("NGEN2_NHV1"); + TwoWindingsTransformer nhv2Nload = network.getTwoWindingsTransformer("NHV2_NLOAD"); + + assertReactivePowerEquals(-250, gen2.getTerminal()); // GEN is correctly limited to 250 MVar + assertReactivePowerEquals(250, ngen2Nhv1.getTerminal1()); + assertReactivePowerEquals(-200, nhv2Nload.getTerminal2()); + + // verify convergence using exact jacobian + knitroParams.setGradientComputationMode(1); + result = loadFlowRunner.run(network, parameters); + assertTrue(result.isFullyConverged()); + + assertReactivePowerEquals(-250, gen2.getTerminal()); // GEN is correctly limited to 250 MVar + assertReactivePowerEquals(250, ngen2Nhv1.getTerminal1()); + assertReactivePowerEquals(-200, nhv2Nload.getTerminal2()); + } + + /** + * Modification of the network to urge a PV-PQ switch and set the new PQ bus to the upper bound on reactive power + */ + @Test + void testReacLimEurostagQup() { + Network network = modifiedEurostagFactory(0, 100); + + // verify convergence using finite differences + LoadFlowResult result = loadFlowRunner.run(network, parameters); + assertTrue(result.isFullyConverged()); + + Generator gen = network.getGenerator("GEN"); + Generator gen2 = network.getGenerator("GEN2"); + TwoWindingsTransformer ngen2Nhv1 = network.getTwoWindingsTransformer("NGEN2_NHV1"); + TwoWindingsTransformer nhv2Nload = network.getTwoWindingsTransformer("NHV2_NLOAD"); + + assertReactivePowerEquals(-280, gen.getTerminal()); + assertReactivePowerEquals(-100, gen2.getTerminal()); // GEN is correctly limited to 100 MVar + assertReactivePowerEquals(100, ngen2Nhv1.getTerminal1()); + assertReactivePowerEquals(-200, nhv2Nload.getTerminal2()); + + // verify convergence using exact jacobian + knitroParams.setGradientComputationMode(1); + result = loadFlowRunner.run(network, parameters); + assertTrue(result.isFullyConverged()); + + assertReactivePowerEquals(-280, gen.getTerminal()); + assertReactivePowerEquals(-100, gen2.getTerminal()); // GEN is correctly limited to 100 MVar + assertReactivePowerEquals(100, ngen2Nhv1.getTerminal1()); + assertReactivePowerEquals(-200, nhv2Nload.getTerminal2()); + } + + /** + * Case of a bus containing a load and a generator. + */ + @Test + void testReacLimEurostagQupWithLoad() { + Network network = modifiedEurostagFactory(0, 100); + + // add reactive load to network + network.getVoltageLevel("VLGEN2") + .newLoad() + .setId("LOAD2") + .setBus("NGEN2") + .setConnectableBus("NGEN2") + .setP0(0.0) + .setQ0(30.0) + .add(); + + // verify convergence using finite differences + LoadFlowResult result = loadFlowRunner.run(network, parameters); + assertTrue(result.isFullyConverged()); + + Generator gen2 = network.getGenerator("GEN2"); + TwoWindingsTransformer ngen2Nhv1 = network.getTwoWindingsTransformer("NGEN2_NHV1"); + TwoWindingsTransformer nhv2Nload = network.getTwoWindingsTransformer("NHV2_NLOAD"); + + assertReactivePowerEquals(-100, gen2.getTerminal()); // GEN is correctly limited to 100 MVar + assertReactivePowerEquals(70, ngen2Nhv1.getTerminal1()); + assertReactivePowerEquals(-200, nhv2Nload.getTerminal2()); + + // verify convergence using exact jacobian + knitroParams.setGradientComputationMode(1); + result = loadFlowRunner.run(network, parameters); + assertTrue(result.isFullyConverged()); + + assertReactivePowerEquals(-100, gen2.getTerminal()); // GEN is correctly limited to 100 MVar + assertReactivePowerEquals(70, ngen2Nhv1.getTerminal1()); + assertReactivePowerEquals(-200, nhv2Nload.getTerminal2()); + } + + /** + * Case of a bus containing two generators. + */ + @Test + void testReacLimEurostagQupWithGen() { + Network network = modifiedEurostagFactory(0, 100); + + // add second generator + Generator gen2Bis = network.getVoltageLevel("VLGEN2") + .newGenerator() + .setId("GEN2BIS") + .setBus("NGEN2") + .setConnectableBus("NGEN2") + .setMinP(-9999.99) + .setMaxP(9999.99) + .setVoltageRegulatorOn(true) + .setTargetV(24.5) + .setTargetP(50) + .add(); + + // add second generator reactive limits + gen2Bis.newMinMaxReactiveLimits() + .setMinQ(0) + .setMaxQ(40) + .add(); + + // verify convergence using finite differences + LoadFlowResult result = loadFlowRunner.run(network, parameters); + assertTrue(result.isFullyConverged()); + + Generator gen = network.getGenerator("GEN"); + Generator gen2 = network.getGenerator("GEN2"); + TwoWindingsTransformer ngen2Nhv1 = network.getTwoWindingsTransformer("NGEN2_NHV1"); + TwoWindingsTransformer nhv2Nload = network.getTwoWindingsTransformer("NHV2_NLOAD"); + + assertReactivePowerEquals(-122.735, gen.getTerminal()); + assertReactivePowerEquals(-100, gen2.getTerminal()); // GEN is correctly limited to 100 MVar + assertReactivePowerEquals(140.0, ngen2Nhv1.getTerminal1()); + assertReactivePowerEquals(-200, nhv2Nload.getTerminal2()); + + // verify convergence using exact jacobian + knitroParams.setGradientComputationMode(1); + result = loadFlowRunner.run(network, parameters); + assertTrue(result.isFullyConverged()); + + assertReactivePowerEquals(-122.735, gen.getTerminal()); + assertReactivePowerEquals(-100, gen2.getTerminal()); // GEN is correctly limited to 100 MVar + assertReactivePowerEquals(140.0, ngen2Nhv1.getTerminal1()); + assertReactivePowerEquals(-200, nhv2Nload.getTerminal2()); + } + + @Test + void testReacLimIeee14() { + Network network = IeeeCdfNetworkFactory.create14(); + LoadFlowResult result = loadFlowRunner.run(network, parameters); + assertTrue(result.isFullyConverged(), "Not Fully Converged"); + } + + @Test + void testReacLimIeee30() { + Network network = IeeeCdfNetworkFactory.create30(); + LoadFlowResult result = loadFlowRunner.run(network, parameters); + assertTrue(result.isFullyConverged(), "Not Fully Converged"); + } + + @Test + void testExactDenseJacobianException() { + // configuration to run generator reactive limits Knitro solver with exact dense Jacobian computation mode + knitroParams.setGradientComputationMode(1) + .setGradientUserRoutine(1); + parameters.addExtension(KnitroLoadFlowParameters.class, knitroParams); + Network network = IeeeCdfNetworkFactory.create14(); + CompletionException e = assertThrows(CompletionException.class, () -> loadFlowRunner.run(network, parameters)); + assertEquals("com.powsybl.commons.PowsyblException: Knitro generator reactive limits is incompatible with exact dense jacobian computation mode: gradientUserRoutine KnitroLoadFlowParameters should be switched to 1", + e.getMessage()); + } + + @Test + void testUseReactiveLimitsOuterLoopException() { + parameters.setUseReactiveLimits(true); + Network network = IeeeCdfNetworkFactory.create14(); + CompletionException e = assertThrows(CompletionException.class, () -> loadFlowRunner.run(network, parameters)); + assertEquals("com.powsybl.commons.PowsyblException: Knitro generator reactive limits and reactive limits outer loop cannot work simultaneously: useReactiveLimits LoadFlowParameter should be switched to false", + e.getMessage()); + } +}