diff --git a/pom.xml b/pom.xml index a34073e..060766c 100644 --- a/pom.xml +++ b/pom.xml @@ -50,7 +50,7 @@ 17 6.6.0 - 1.14.1 + 1.14.1-SNAPSHOTYA 14.2.0 0.7.0 diff --git a/src/main/java/com/powsybl/openloadflow/knitro/solver/KnitroLoadFlowParameters.java b/src/main/java/com/powsybl/openloadflow/knitro/solver/KnitroLoadFlowParameters.java index e7fc73b..bd6a8d9 100644 --- a/src/main/java/com/powsybl/openloadflow/knitro/solver/KnitroLoadFlowParameters.java +++ b/src/main/java/com/powsybl/openloadflow/knitro/solver/KnitroLoadFlowParameters.java @@ -27,6 +27,15 @@ public class KnitroLoadFlowParameters extends AbstractExtension new ResilientKnitroSolver(network, knitroSolverParameters, equationSystem, j, targetVector, equationVector, parameters.isDetailedReport()); + case REACTIVLIMITS -> + new KnitroSolverReacLim(network, knitroSolverParameters, equationSystem, j, targetVector, equationVector, parameters.isDetailedReport()); }; } 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 e3053e6..128c694 100644 --- a/src/main/java/com/powsybl/openloadflow/knitro/solver/KnitroSolverParameters.java +++ b/src/main/java/com/powsybl/openloadflow/knitro/solver/KnitroSolverParameters.java @@ -22,6 +22,7 @@ */ public class KnitroSolverParameters implements AcSolverParameters { + public static final String DEFAULT_LOG_FILE_NAME = "Logs.txt"; public static final int DEFAULT_GRADIENT_COMPUTATION_MODE = 1; // Specifies how the Jacobian matrix is computed public static final int DEFAULT_GRADIENT_USER_ROUTINE = 2; // If the user chooses to pass the exact Jacobian to knitro, specifies the sparsity pattern for the Jacobian matrix. public static final int DEFAULT_HESSIAN_COMPUTATION_MODE = 6; // Specifies how the Hessian matrix is computed. 6 means that the Hessian is approximated using the L-BFGS method, which is a quasi-Newton method. @@ -31,7 +32,7 @@ public class KnitroSolverParameters implements AcSolverParameters { public static final double DEFAULT_STOPPING_CRITERIA = Math.pow(10, -6); public static final StateVectorScalingMode DEFAULT_STATE_VECTOR_SCALING_MODE = StateVectorScalingMode.NONE; public static final boolean ALWAYS_UPDATE_NETWORK_DEFAULT_VALUE = true; - public static final boolean CHECK_LOAD_FLOW_SOLUTION_DEFAULT_VALUE = true; // If the solver converges, launches a load flow taking into account slack variable results to check if the solution is a real load flow solution + public static final boolean CHECK_LOAD_FLOW_SOLUTION_DEFAULT_VALUE = false; // If the solver converges, launches a load flow taking into account slack variable results to check if the solution is a real load flow solution public static final KnitroSolverType DEFAULT_KNITRO_SOLVER_TYPE = KnitroSolverType.STANDARD; public KnitroSolverType knitroSolverType = DEFAULT_KNITRO_SOLVER_TYPE; private StateVectorScalingMode stateVectorScalingMode = DEFAULT_STATE_VECTOR_SCALING_MODE; @@ -48,6 +49,16 @@ public class KnitroSolverParameters implements AcSolverParameters { private int maxIterations = DEFAULT_MAX_ITERATIONS; private double convEps = DEFAULT_STOPPING_CRITERIA; private int hessianComputationMode = DEFAULT_HESSIAN_COMPUTATION_MODE; + private String logFileName = DEFAULT_LOG_FILE_NAME; + private KnitroWritter knitroWritter; + + public KnitroSolverParameters(KnitroWritter knitroWritter) { + this.knitroWritter = knitroWritter; + } + + public KnitroSolverParameters() { + this.knitroWritter = new KnitroWritter("Logs.txt"); + } public int getGradientComputationMode() { return gradientComputationMode; @@ -202,6 +213,15 @@ public void setKnitroSolverType(KnitroSolverType knitroSolverType) { this.knitroSolverType = Objects.requireNonNull(knitroSolverType); } + public KnitroWritter getKnitroWritter() { + return knitroWritter; + } + + public KnitroSolverParameters setKnitroWritter(KnitroWritter knitroWritter) { + this.knitroWritter = knitroWritter; + return this; + } + @Override public String toString() { return "KnitroSolverParameters(" + @@ -218,6 +238,7 @@ public String toString() { public enum KnitroSolverType { STANDARD, - RESILIENT + RESILIENT, + REACTIVLIMITS, } } diff --git a/src/main/java/com/powsybl/openloadflow/knitro/solver/KnitroSolverReacLim.java b/src/main/java/com/powsybl/openloadflow/knitro/solver/KnitroSolverReacLim.java new file mode 100644 index 0000000..0efeabb --- /dev/null +++ b/src/main/java/com/powsybl/openloadflow/knitro/solver/KnitroSolverReacLim.java @@ -0,0 +1,1461 @@ +/** + * 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.KNEvalFCCallback; +import com.artelys.knitro.api.callbacks.KNEvalGACallback; +import com.artelys.knitro.api.nativelibrary.KNLibrary; +import com.powsybl.commons.PowsyblException; +import com.powsybl.commons.report.ReportNode; +import com.powsybl.math.matrix.SparseMatrix; +import com.powsybl.openloadflow.ac.equations.AcEquationType; +import com.powsybl.openloadflow.ac.equations.AcVariableType; +import com.powsybl.openloadflow.ac.solver.AbstractAcSolver; +import com.powsybl.openloadflow.ac.solver.AcSolverResult; +import com.powsybl.openloadflow.ac.solver.AcSolverStatus; +import com.powsybl.openloadflow.ac.solver.AcSolverUtil; +import com.powsybl.openloadflow.equations.*; +import com.powsybl.openloadflow.network.*; +import com.powsybl.openloadflow.network.util.VoltageInitializer; +import org.apache.commons.lang3.Range; +import org.slf4j.Logger; +import org.slf4j.LoggerFactory; + +import java.util.*; +import java.util.stream.Collectors; +import java.util.stream.IntStream; +import java.util.stream.Stream; + +import static com.artelys.knitro.api.nativelibrary.KNLibrary.KN_get_solve_time_cpu; +import static com.google.common.primitives.Doubles.toArray; +import static com.powsybl.openloadflow.ac.equations.AcEquationType.*; + +/** + * @author Martin Debouté {@literal } + * @author Pierre Arvy {@literal } + * @author Yoann Anezin {@literal } + */ +public class KnitroSolverReacLim extends AbstractAcSolver { + + private static final Logger LOGGER = LoggerFactory.getLogger(KnitroSolverReacLim.class); + private boolean firstIter = true; + private final KnitroWritter knitroWritter; + + // Penalty weights in the objective function + private final double wK = 1.0; + private final double wP = 1.0; + private final double wQ = 1.0; + private final double wV = 10.0; + + // Lambda + private final double lambda = 2.0; + private final double mu = 1.0; + + // Number of Load Flows (LF) variables in the system + private final int numLFVariables; + + // Number of variables including slack variables + private final int numLFandSlackVariables; + // Total number of variables including complementarity constraints' ones + private final int numTotalVariables; + + // 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; + private final int complConstVariables; + + // Starting indices for slack variables in the variable vector + private final int slackStartIndex; + private final int slackPStartIndex; + private final int slackQStartIndex; + private final int slackVStartIndex; + private final int compVarStartIndex; + + // Mappings from global equation indices to local indices by equation type + private final Map pEquationLocalIds; + private final Map qEquationLocalIds; + private final Map vEquationLocalIds; + private static final Map> INDEQUNACTIVEQ = new LinkedHashMap<>(); + + // Mapping of slacked bus + private final ArrayList slackContributions = new ArrayList<>(); + + // Unactivated Equations on reactiv power to deal with + private final List> equationsQBusV; + private final List listElementNumWithQEqUnactivated; + private final Map vSuppEquationLocalIds; + protected KnitroSolverParameters knitroParameters; + + // TIme spent in class + private long time = 0; + + public KnitroSolverReacLim( + LfNetwork network, + KnitroSolverParameters knitroParameters, + EquationSystem equationSystem, + JacobianMatrix jacobian, + TargetVector targetVector, + EquationVector equationVector, + boolean detailedReport) { + super(network, equationSystem, jacobian, targetVector, equationVector, detailedReport); + long start = System.nanoTime(); + this.knitroParameters = knitroParameters; + this.knitroWritter = knitroParameters.getKnitroWritter(); + + this.numLFVariables = equationSystem.getIndex().getSortedVariablesToFind().size(); + + List> sortedEquations = equationSystem.getIndex().getSortedEquationsToSolve(); + + // Count number of classic LF 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.numLFandSlackVariables = numLFVariables + numSlackVariables; + + this.slackStartIndex = numLFVariables; + this.slackPStartIndex = slackStartIndex; + this.slackQStartIndex = slackPStartIndex + 2 * numPEquations; + this.slackVStartIndex = slackQStartIndex + 2 * numQEquations; + this.compVarStartIndex = slackVStartIndex + 2 * numVEquations; + + // The new equation system implemented here duplicates and modifies V equations + // It also adds two equations on Q on each PV bus + // First we need to collect those Q equations + List> activeEquationsV = sortedEquations.stream() + .filter(e -> e.getType() == AcEquationType.BUS_TARGET_V).toList(); + List listBusesWithVEq = activeEquationsV.stream() + .map(e -> e.getTerms().get(0).getElementNum()).toList(); + List> equationQBusElementNum = new ArrayList<>(); + for (int elementNum : listBusesWithVEq) { + List> listEqElementNum = equationSystem.getEquations(ElementType.BUS, elementNum); + equationQBusElementNum.addAll(listEqElementNum.stream().filter(e -> + e.getType() == BUS_TARGET_Q).toList()); + } + + // Are taking into account only buses with limits on reactive power, the others are left as in the initial model + List listBusesWithQEqToAdd = new ArrayList<>(); + List> equationQBusVToAdd = new ArrayList<>(); + for (Equation equation : equationQBusElementNum) { + boolean limitedOnQ = true; + for (LfGenerator g : network.getBuses().get(equation.getElement(network).get().getNum()).getGenerators()) { + if (g.getMaxQ() >= 1.7976931348623156E306 || g.getMinQ() <= -1.7976931348623156E306) { + limitedOnQ = false; + } + } + if (limitedOnQ) { + equationQBusVToAdd.add(equation); + listBusesWithQEqToAdd.add(equation.getElementNum()); + } + } + + // 3 new variables are used on both V equations modified, and 2 on the Q equations listed above + this.complConstVariables = equationQBusVToAdd.size() * 5; + + this.numTotalVariables = numLFandSlackVariables + complConstVariables; + this.equationsQBusV = Stream.concat(equationQBusVToAdd.stream(), + equationQBusVToAdd.stream()).toList(); //Duplication to get b_low and b_up eq + this.listElementNumWithQEqUnactivated = listBusesWithQEqToAdd; + + // Map equations to local indices + this.pEquationLocalIds = new HashMap<>(); + this.qEquationLocalIds = new HashMap<>(); + this.vEquationLocalIds = new HashMap<>(); + this.vSuppEquationLocalIds = new HashMap<>(); + + int pCounter = 0; + int qCounter = 0; + int vCounter = 0; + int vSuppCounter = 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++); + break; + case BUS_TARGET_Q: + qEquationLocalIds.put(i, qCounter++); + break; + case BUS_TARGET_V: + vEquationLocalIds.put(i, vCounter++); + if (listElementNumWithQEqUnactivated.contains(sortedEquations.get(i).getElementNum())) { + vSuppEquationLocalIds.put(i, vSuppCounter++); + } + break; + } + } + knitroWritter.write("Poids de la fonction objectif : wK = " + wK + ", wP = " + wP + ", wQ = " + wQ + ", wV =" + wV, true); + knitroWritter.write("Nombre de Variables de LoadFLow : " + numLFVariables, true); + knitroWritter.write("Nombre de Variables de Slacks : " + 2 * (numPEquations + numQEquations + numVEquations), true); + knitroWritter.write("Nombre de Variables de Complémentarités Initialement Prévues : " + 5 * numVEquations, true); + knitroWritter.write("Nombre total de Variables : " + numTotalVariables, true); + + long end = System.nanoTime(); + time += end - start; + } + + /** + * Returns the name of the solver. + */ + @Override + public String getName() { + return "Knitro Reactive Limits Solver"; + } + + /** + * Logs the Knitro status and its corresponding AcSolverStatus. + * + * @param status the KnitroStatus to log + */ + private void logKnitroStatus(KnitroStatus status) { + LOGGER.info("Knitro Status: {}", status); + } + + /** + * Builds the row and column index lists corresponding to the dense Jacobian structure, + * assuming each non-linear constraint is derived with respect to every variable. + * + * @param numVars Total number of variables. + * @param listNonLinearConsts List of non-linear constraint indices. + * @param listNonZerosCtsDense Output list to receive constraint indices for non-zero Jacobian entries. + * @param listNonZerosVarsDense Output list to receive variable indices for non-zero Jacobian entries. + */ + public void buildDenseJacobianMatrix( + int numVars, + List listNonLinearConsts, + List listNonZerosCtsDense, + List listNonZerosVarsDense) { + long start = System.nanoTime(); + // Each non-linear constraint will have a partial derivative with respect to every variable + for (Integer constraintId : listNonLinearConsts) { + for (int varIndex = 0; varIndex < numVars; varIndex++) { + listNonZerosCtsDense.add(constraintId); + } + } + + // We repeat the full list of variables for each non-linear constraint + List variableIndices = IntStream.range(0, numVars).boxed().toList(); + for (int i = 0; i < listNonLinearConsts.size(); i++) { + listNonZerosVarsDense.addAll(variableIndices); + } + long end = System.nanoTime(); + time += end - start; + } + + /** + * Builds the sparse Jacobian matrix by identifying non-zero entries + * for each non-linear constraint, including contributions from slack variables. + * + * @param sortedEquationsToSolve Ordered list of equations to solve. + * @param nonLinearConstraintIds Indices of non-linear constraints within the sorted equation list. + * @param jacobianRowIndices Output: row indices (constraints) of non-zero Jacobian entries. + * @param jacobianColumnIndices Output: column indices (variables) of non-zero Jacobian entries. + */ + public void buildSparseJacobianMatrix( + List> sortedEquationsToSolve, + List nonLinearConstraintIds, + List jacobianRowIndices, + List jacobianColumnIndices) { + long start = System.nanoTime(); + int numberLFEq = sortedEquationsToSolve.size() - 3 * vSuppEquationLocalIds.size(); + + for (Integer constraintIndex : nonLinearConstraintIds) { + Equation equation = sortedEquationsToSolve.get(constraintIndex); + AcEquationType equationType = equation.getType(); + + // Collect all variable indices involved in the current equation + Set involvedVariables = equation.getTerms().stream() + .flatMap(term -> term.getVariables().stream()) + .map(Variable::getRow) + .collect(Collectors.toCollection(TreeSet::new)); + + // Add slack variables if the constraint type has them + 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) { + 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); + }; + involvedVariables.add(slackBaseIndex + 2 * slackStart); // Sm + involvedVariables.add(slackBaseIndex + 2 * slackStart + 1); // Sp + } + + // Add complementarity constraints' variables if the constraint type has them + int compVarStart; + // Case of inactive Q equations + if (equationType == BUS_TARGET_Q && !equation.isActive()) { + compVarStart = vSuppEquationLocalIds.getOrDefault(equationSystem.getIndex().getSortedEquationsToSolve() + .indexOf(equationSystem.getEquations(ElementType.BUS, equation.getElementNum()).stream() + .filter(e -> e.getType() == BUS_TARGET_V).toList() + .get(0)), -1); + if (constraintIndex < numberLFEq + 2 * vSuppEquationLocalIds.size()) { + involvedVariables.add(compVarStartIndex + 5 * compVarStart + 2); // b_low + } else { + involvedVariables.add(compVarStartIndex + 5 * compVarStart + 3); // b_up + } + } + + // Add one entry for each non-zero (constraintIndex, variableIndex) + jacobianColumnIndices.addAll(involvedVariables); + jacobianRowIndices.addAll(Collections.nCopies(involvedVariables.size(), constraintIndex)); + + long end = System.nanoTime(); + time += end - start; + } + } + + /** + * Sets Knitro solver parameters based on the provided KnitroSolverParameters object. + * + * @param solver The Knitro solver instance to configure. + * @throws KNException if Knitro fails to accept a parameter. + */ + private void setSolverParameters(KNSolver solver) throws KNException { + long start = System.nanoTime(); + LOGGER.info("Configuring Knitro solver parameters..."); + + solver.setParam(KNConstants.KN_PARAM_GRADOPT, knitroParameters.getGradientComputationMode()); + solver.setParam(KNConstants.KN_PARAM_FEASTOL, knitroParameters.getConvEps()); + solver.setParam(KNConstants.KN_PARAM_MAXIT, knitroParameters.getMaxIterations()); + solver.setParam(KNConstants.KN_PARAM_HESSOPT, knitroParameters.getHessianComputationMode()); +// solver.setParam(KNConstants.KN_PARAM_SOLTYPE, KNConstants.KN_SOLTYPE_BESTFEAS); +// solver.setParam(KNConstants.KN_PARAM_OUTMODE, KNConstants.KN_OUTMODE_BOTH); + solver.setParam(KNConstants.KN_PARAM_OPTTOL, 1.0e-3); + solver.setParam(KNConstants.KN_PARAM_OPTTOLABS, 1.0e-1); + solver.setParam(KNConstants.KN_PARAM_OUTLEV, 3); +// solver.setParam(KNConstants.KN_PARAM_NUMTHREADS, 1); +// solver.setParam(KNConstants.KN_PARAM_BAR_MPEC_HEURISTIC, 1); + + LOGGER.info("Knitro parameters set: GRADOPT={}, HESSOPT={}, FEASTOL={}, MAXIT={}", + knitroParameters.getGradientComputationMode(), + knitroParameters.getHessianComputationMode(), + knitroParameters.getConvEps(), + knitroParameters.getMaxIterations()); + long end = System.nanoTime(); + time += end - start; + } + + @Override + public AcSolverResult run(VoltageInitializer voltageInitializer, ReportNode reportNode) { + long start = System.nanoTime(); + int nbIterations; + AcSolverStatus solverStatus; + ResilientReacLimKnitroProblem problemInstance; + + try { + problemInstance = new ResilientReacLimKnitroProblem(network, equationSystem, targetVector, j, voltageInitializer); + } catch (KNException e) { + throw new PowsyblException("Exception while building Knitro problem", e); + } + + try { + KNSolver solver = new KNSolver(problemInstance); + solver.initProblem(); + setSolverParameters(solver); + solver.solve(); + KNSolution solution = solver.getSolution(); + List constraintValues = solver.getConstraintValues(); + List x = solution.getX(); + List lambda2 = solution.getLambda(); + + solverStatus = KnitroStatus.fromStatusCode(solution.getStatus()).toAcSolverStatus(); + logKnitroStatus(KnitroStatus.fromStatusCode(solution.getStatus())); + nbIterations = solver.getNumberIters(); + + LOGGER.info("==== Solution Summary ===="); + LOGGER.info("Objective value = {}", solution.getObjValue()); + LOGGER.info("Feasibility violation = {}", solution.getFeasError()); + LOGGER.info("Optimality violation = {}", solver.getAbsOptError()); + knitroWritter.write("==== Solution Summary ====", true); + knitroWritter.write("Objective value = " + solution.getObjValue(), true); + knitroWritter.write("Feasibility violation = " + solution.getFeasError(), true); + knitroWritter.write("Optimality violation = " + solver.getAbsOptError(), true); + + // Log primal solution + LOGGER.debug("==== Optimal variables ===="); + for (int i = 0; i < x.size(); i++) { + LOGGER.debug(" x[{}] = {}", i, x.get(i)); + } + + LOGGER.debug("==== Constraint values ===="); + for (int i = 0; i < problemInstance.getNumCons(); i++) { + LOGGER.debug(" c[{}] = {} (λ = {})", i, constraintValues.get(i), lambda2.get(i)); + } + + LOGGER.debug("==== Constraint violations ===="); + for (int i = 0; i < problemInstance.getNumCons(); i++) { + LOGGER.debug(" violation[{}] = {}", i, solver.getConViol(i)); + } + + // ========== 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, wK * wP, lambda, mu); + double penaltyQ = computeSlackPenalty(x, slackQStartIndex, numQEquations, wK * wQ, lambda, mu); + double penaltyV = computeSlackPenalty(x, slackVStartIndex, numVEquations, wV, lambda, mu); + 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); + + knitroWritter.write("Penalty P = " + penaltyP, true); + knitroWritter.write("Penalty Q = " + penaltyQ, true); + knitroWritter.write("Penalty V = " + penaltyV, true); + knitroWritter.write("Total penalty = " + totalPenalty, true); + + LOGGER.info("=== Switches Done==="); + checkSwitchesDone(x, compVarStartIndex, complConstVariables / 5); + + // ========== Network Update ========== + // Update the network values if the solver converged or if the network should always be updated + if (solverStatus == AcSolverStatus.CONVERGED || knitroParameters.isAlwaysUpdateNetwork()) { + //Update the state vector with the solution + equationSystem.getStateVector().set(toArray(x)); + for (Equation equation : equationSystem.getEquations()) { + for (EquationTerm term : equation.getTerms()) { + term.setStateVector(equationSystem.getStateVector()); + } + } + AcSolverUtil.updateNetwork(network, equationSystem); + } + + } catch (KNException e) { + throw new PowsyblException("Exception while solving with Knitro", e); + } + + double slackBusMismatch = network.getSlackBuses().stream() + .mapToDouble(LfBus::getMismatchP) + .sum(); + + long end = System.nanoTime(); + time += end - start; +// knitroWritter.write("Temps passé dans la classe KnitroSolverReacLim = " + time * 1e-9, true); + return new AcSolverResult(solverStatus, nbIterations, slackBusMismatch); + } + + private void logSlackValues(String type, int startIndex, int count, List x) { + long start = System.nanoTime(); + KnitroWritter slackPWritter = new KnitroWritter("D:\\Documents\\Slacks\\SlacksP.txt"); + KnitroWritter slackQWritter = new KnitroWritter("D:\\Documents\\Slacks\\SlacksQ.txt"); + KnitroWritter slackVWritter = new KnitroWritter("D:\\Documents\\Slacks\\SlacksV.txt"); + final double threshold = 1e-6; // Threshold for significant slack values + final double sbase = 100.0; // Base power in MVA + + LOGGER.info("==== Slack diagnostics for {} (p.u. and physical units) ====", type); + boolean firstIterP = true; + boolean firstIterQ = true; + boolean firstIterV = true; + 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; + + if (Math.abs(epsilon) <= threshold) { + continue; + } + + String name = getSlackVariableBusName(i, type); + String interpretation; + + 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); + continue; + } + interpretation = String.format("ΔV = %.4f p.u. (%.1f kV)", epsilon, epsilon * bus.getNominalV()); + } + default -> interpretation = "Unknown slack type"; + } + + slackContributions.add(new ResilientKnitroSolver.SlackKey(type, name, epsilon)); + String msg = String.format("Slack %s[ %s ] → Sm = %.4f, Sp = %.4f → %s", type, name, sm, sp, interpretation); + LOGGER.info(msg); + knitroWritter.write(msg, true); + switch (type) { + case "P": + slackPWritter.write(name, !firstIterP); + slackPWritter.write(String.format("%.4f", epsilon), true); + firstIterP = false; + break; + case "Q": + slackQWritter.write(name, !firstIterQ); + slackQWritter.write(String.format("%.4f", epsilon), true); + firstIterQ = false; + break; + case "V": + var bus = network.getBusById(name); + slackVWritter.write(name, !firstIterV); + slackVWritter.write(String.format("%.4f", epsilon * bus.getNominalV()), true); + firstIterV = false; + break; + } + } + long end = System.nanoTime(); + time += end - start; + } + + private String getSlackVariableBusName(Integer index, String type) { + long start = System.nanoTime(); + 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 RuntimeException("Variable index associated with slack variable " + type + "was not found"); + } + + LfBus bus = network.getBus(equationSystem.getIndex().getSortedEquationsToSolve().get(varIndex).getElementNum()); + + long end = System.nanoTime(); + time += end - start; + return bus.getId(); + } + + private double computeSlackPenalty(List x, int startIndex, int count, double weight, double lambda, double mu) { + long start = System.nanoTime(); + 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 * mu * (diff * diff); // Quadratic terms + penalty += weight * lambda * (sp + sm); // Linear terms + } + long end = System.nanoTime(); + time += end - start; + return penalty; + } + + /** + * Inform all switches PV -> PQ done in the solution found + * @param x current network's estate + * @param startIndex first index of complementarity constraints variables in x + * @param count number of b_low / b_up different variables + */ + private void checkSwitchesDone(List x, int startIndex, int count) { + long start = System.nanoTime(); + int nombreSwitches = 0; + 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); + String bus = equationSystem.getIndex() + .getSortedEquationsToSolve().stream().filter(e -> + e.getType() == BUS_TARGET_V).toList().get(i).getElement(network).get().getId(); + if (Math.abs(bLow) < 1E-4 && !(vInf < 1E-4 && vSup < 1E-4)) { + nombreSwitches++; + LOGGER.info("Switch PV -> PQ on bus {}, Q set at Qmin", bus); + knitroWritter.write("Switch PV -> PQ on bus " + bus + ", Q set at Qmin", true); + + } else if (Math.abs(bUp) < 1E-4 && !(vInf < 1E-4 && vSup < 1E-4)) { + nombreSwitches++; + LOGGER.info("Switch PV -> PQ on bus {}, Q set at Qmax", bus); + knitroWritter.write("Switch PV -> PQ on bus " + bus + ", Q set at Qmax", true); + + } + } + knitroWritter.write("Nombre total de switches : " + nombreSwitches, true); + long end = System.nanoTime(); + time += end - start; + } + + /** + * Returns the sparsity pattern of the hessian matrix associated with the problem. + * + * @param nonlinearConstraintIndexes A list of the indexes of non-linear equations. + * @return row and column coordinates of non-zero entries in the hessian matrix. + */ + private AbstractMap.SimpleEntry, List> getHessNnzRowsAndCols(List nonlinearConstraintIndexes, List> equationsToSolve) { + record NnzCoordinates(int iRow, int iCol) { + } + long start = System.nanoTime(); + Set hessianEntries = new LinkedHashSet<>(); + + // Non-linear constraints contributions in the hessian matrix + for (int index : nonlinearConstraintIndexes) { + if (index < equationsToSolve.size()) { + Equation equation = equationsToSolve.get(index); + for (EquationTerm term : equation.getTerms()) { + for (Variable var1 : term.getVariables()) { + int i = var1.getRow(); + for (Variable var2 : term.getVariables()) { + int j = var2.getRow(); + if (j >= i) { + hessianEntries.add(new NnzCoordinates(i, j)); + } + } + } + } + } + } + + // Slacks variables contributions in the objective function + for (int iSlack = slackStartIndex; iSlack < compVarStartIndex; iSlack++) { + hessianEntries.add(new NnzCoordinates(iSlack, iSlack)); + if (((iSlack - slackStartIndex) & 1) == 0) { + hessianEntries.add(new NnzCoordinates(iSlack, iSlack + 1)); + } + } + + // Sort the entries by row and column indices + hessianEntries = hessianEntries.stream() + .sorted(Comparator.comparingInt(NnzCoordinates::iRow).thenComparingInt(NnzCoordinates::iCol)) + .collect(Collectors.toCollection(LinkedHashSet::new)); + + List hessRows = new ArrayList<>(); + List hessCols = new ArrayList<>(); + for (NnzCoordinates entry : hessianEntries) { + hessRows.add(entry.iRow()); + hessCols.add(entry.iCol()); + } + + long end = System.nanoTime(); + time += end - start; + return new AbstractMap.SimpleEntry<>(hessRows, hessCols); + } + + /** + * Enum representing specific status codes returned by the Knitro solver, + * grouped either individually or by ranges, and mapped to corresponding {@link AcSolverStatus} values. + * This mapping allows a more fine-grained interpretation of solver termination reasons, + * distinguishing between convergence, infeasibility, modeling errors, evaluation issues, etc... + */ + public enum KnitroStatus { + + /** + * Successful convergence to a local optimum. + */ + CONVERGED_TO_LOCAL_OPTIMUM(0, 0, AcSolverStatus.CONVERGED), + + /** + * Converged to a feasible but not necessarily optimal solution. + */ + CONVERGED_TO_FEASIBLE_APPROXIMATE_SOLUTION(-199, -100, AcSolverStatus.CONVERGED), + + /** + * Solver terminated at an infeasible point. + */ + TERMINATED_AT_INFEASIBLE_POINT(-299, -200, AcSolverStatus.SOLVER_FAILED), + + /** + * The problem was detected as unbounded. + */ + PROBLEM_UNBOUNDED(-399, -300, AcSolverStatus.SOLVER_FAILED), + + /** + * Optimization stopped due to reaching iteration or time limits. + */ + TERMINATED_DUE_TO_PRE_DEFINED_LIMIT(-499, -400, AcSolverStatus.MAX_ITERATION_REACHED), + + /** + * Failure in a user-defined callback function. + */ + CALLBACK_ERROR(-500, -500, AcSolverStatus.SOLVER_FAILED), + + /** + * Internal LP solver failure in active-set method. + */ + LP_SOLVER_ERROR(-501, -501, AcSolverStatus.SOLVER_FAILED), + + /** + * Evaluation failure (e.g., division by zero or invalid sqrt). + */ + EVALUATION_ERROR(-502, -502, AcSolverStatus.SOLVER_FAILED), + + /** + * Insufficient memory to solve the problem. + */ + OUT_OF_MEMORY(-503, -503, AcSolverStatus.SOLVER_FAILED), + + /** + * Solver was stopped manually by the user. + */ + USER_TERMINATION(-504, -504, AcSolverStatus.SOLVER_FAILED), + + /** + * File open error when trying to read input. + */ + INPUT_FILE_ERROR(-505, -505, AcSolverStatus.SOLVER_FAILED), + + /** + * Modeling error: invalid variable/constraint setup. + */ + MODEL_DEFINITION_ERROR(-530, -506, AcSolverStatus.SOLVER_FAILED), + + /** + * Internal Knitro error – contact support. + */ + INTERNAL_ERROR(-600, -600, AcSolverStatus.SOLVER_FAILED), + + /** + * Fallback for unknown status codes. + */ + UNKNOWN_STATUS(Integer.MIN_VALUE, Integer.MAX_VALUE, AcSolverStatus.SOLVER_FAILED); + + private final Range codeRange; + private final AcSolverStatus mappedStatus; + + /** + * Constructs a KnitroStatus with a range of status codes and the associated AcSolverStatus. + * + * @param min The minimum status code value (inclusive). + * @param max The maximum status code value (inclusive). + * @param mappedStatus The corresponding AcSolverStatus. + */ + KnitroStatus(int min, int max, AcSolverStatus mappedStatus) { + this.codeRange = Range.of(min, max); + this.mappedStatus = mappedStatus; + } + + /** + * Returns the KnitroStatus corresponding to the given status code. + * + * @param statusCode the status code returned by Knitro + * @return the matching KnitroStatus enum constant, or UNKNOWN_STATUS if unknown + */ + public static KnitroStatus fromStatusCode(int statusCode) { + return Arrays.stream(KnitroStatus.values()) + .filter(status -> status.codeRange.contains(statusCode)) + .findFirst() + .orElse(UNKNOWN_STATUS); + } + + /** + * Returns the {@link AcSolverStatus} associated with this KnitroStatus. + * + * @return the corresponding AcSolverStatus + */ + public AcSolverStatus toAcSolverStatus() { + return mappedStatus; + } + } + + record SlackKey(String type, String busId, double contribution) { + } + + private final class ResilientReacLimKnitroProblem extends KNProblem { + + /** + * Knitro problem definition including: + * - Initialization of variables (types, bounds, initial state) + * - Definition of linear and non-linear constraints + * - Objective function setup + * - Jacobian matrix setup for Knitro + */ + private ResilientReacLimKnitroProblem( + LfNetwork network, + EquationSystem equationSystem, + TargetVector targetVector, + JacobianMatrix jacobianMatrix, + VoltageInitializer voltageInitializer) throws KNException { + // =============== Variable Initialization =============== + super(numTotalVariables, equationSystem.getIndex().getSortedEquationsToSolve().size() + + 3 * complConstVariables / 5); + long start = System.nanoTime(); + + // Variable types (all continuous), bounds, and initial values + List 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)); + List scalingFactors = new ArrayList<>(Collections.nCopies(numTotalVariables, 1.0)); + List scalingCenters = new ArrayList<>(Collections.nCopies(numTotalVariables, 0.0)); + List initialValues = new ArrayList<>(Collections.nCopies(numTotalVariables, 0.0)); + + setVarTypes(variableTypes); + + // Compute initial voltage state using the given initializer + AcSolverUtil.initStateVector(network, equationSystem, voltageInitializer); + for (int i = 0; i < numLFVariables; i++) { + initialValues.set(i, equationSystem.getStateVector().get(i)); + } + + // Initialize slack variables (≥ 0, initial value = 0), scale P and Q slacks + boolean scaled = false; + for (int i = slackStartIndex; i < numLFandSlackVariables; i++) { + lowerBounds.set(i, 0.0); +// upperBounds.set(i, 0.0); +// if (i >= slackQStartIndex) { +// upperBounds.set(i, 0.0); +// } + + if (i < slackVStartIndex) { + scalingFactors.set(i, 1e-2); + scaled = true; + } + + if (i >= slackVStartIndex) { +// upperBounds.set(i,0.0); + } + if (scaled && firstIter) { + knitroWritter.write("Scaling value sur les slacks P et Q : " + 1e-2, true); + firstIter = false; + } else if (firstIter) { + knitroWritter.write("No Scaling applied", true); + firstIter = false; + } + } + + // Set bounds for voltage variables based on Knitro parameters + for (int i = 0; i < numLFVariables; i++) { + if (equationSystem.getIndex().getSortedVariablesToFind().get(i).getType() == AcVariableType.BUS_V) { + lowerBounds.set(i, knitroParameters.getLowerVoltageBound()); + upperBounds.set(i, knitroParameters.getUpperVoltageBound()); + } + } + + // Set bounds for complementarity variables (≥ 0) + for (int i = 0; i < complConstVariables / 5; 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); + } + + LOGGER.info("Voltage initialization strategy: {}", voltageInitializer); + + int n = numTotalVariables; + ArrayList list = new ArrayList<>(n); + for (int i = 0; i < n; i++) { + list.add(i); + } + + // Set bounds and initial state + setVarLoBnds(lowerBounds); + setVarUpBnds(upperBounds); + setVarScaleFactors(new KNSparseVector<>(list, scalingFactors)); + setVarScaleCenters(new KNSparseVector<>(list, scalingCenters)); + setXInitial(initialValues); + LOGGER.info("Variables initialization complete!"); + + List> activeConstraints = equationSystem.getIndex().getSortedEquationsToSolve(); + + // Linear and nonlinear constraints (the latter are deferred to callback) + NonLinearExternalSolverUtils solverUtils = new NonLinearExternalSolverUtils(); + + List nonlinearConstraintIndexes = new ArrayList<>(); // contains the indexes of all non-linear constraints + List> completeEquationsToSolve = new ArrayList<>(activeConstraints); + List wholeTargetVector = new ArrayList<>(Arrays.stream(targetVector.getArray()).boxed().toList()); + for (int equationId = 0; equationId < activeConstraints.size(); equationId++) { + addActivatedConstraints(equationId, activeConstraints, solverUtils, nonlinearConstraintIndexes, + completeEquationsToSolve, targetVector, wholeTargetVector, listElementNumWithQEqUnactivated); // Add Linear constraints, index nonLinear ones and get target values + } + int totalActiveConstraints = completeEquationsToSolve.size(); + completeEquationsToSolve.addAll(equationsQBusV); + + // Set Target Q on the unactive equations added + for (int equationId = totalActiveConstraints; equationId < completeEquationsToSolve.size(); equationId++) { + Equation equation = completeEquationsToSolve.get(equationId); + LfBus b = network.getBus(equation.getElementNum()); + if (equationId - totalActiveConstraints < equationsQBusV.size() / 2) { + wholeTargetVector.add(b.getMinQ() - b.getLoadTargetQ()); + } else { + wholeTargetVector.add(b.getMaxQ() - b.getLoadTargetQ()); + } + nonlinearConstraintIndexes.add(equationId); + INDEQUNACTIVEQ.put(equationId, equation); + } + + int numConstraints = completeEquationsToSolve.size(); + LOGGER.info("Defined {} constraints", numConstraints); + setMainCallbackCstIndexes(nonlinearConstraintIndexes); + setConEqBnds(wholeTargetVector); + + // =============== Declaration of Complementarity Constraints =============== + List listTypeVar = new ArrayList<>(Collections.nCopies(2 * complConstVariables / 5, KNConstants.KN_CCTYPE_VARVAR)); + List bVarList = new ArrayList<>(); // b_up, b_low + List vInfSuppList = new ArrayList<>(); // V_inf, V_sup + + for (int i = 0; i < complConstVariables / 5; 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 + } + + setCompConstraintsTypes(listTypeVar); + setCompConstraintsParts(bVarList, vInfSuppList); // b_low compl with V_sup and b_up compl with V_inf + + // =============== Objective Function =============== + List quadRows = new ArrayList<>(); + List quadCols = new ArrayList<>(); + List quadCoefs = new ArrayList<>(); + + List linIndexes = new ArrayList<>(); + List linCoefs = new ArrayList<>(); + + // Slack penalty terms: (Sp - Sm)^2 = Sp^2 + Sm^2 - 2*Sp*Sm + linear terms from the absolute value + addSlackObjectiveTerms(numPEquations, slackPStartIndex, wK * wP, lambda, mu, quadRows, quadCols, quadCoefs, linIndexes, linCoefs); + addSlackObjectiveTerms(numQEquations, slackQStartIndex, wK * wQ, lambda, mu, quadRows, quadCols, quadCoefs, linIndexes, linCoefs); + addSlackObjectiveTerms(numVEquations, slackVStartIndex, wV, lambda, mu, quadRows, quadCols, quadCoefs, linIndexes, linCoefs); + + setObjectiveQuadraticPart(quadRows, quadCols, quadCoefs); + setObjectiveLinearPart(linIndexes, linCoefs); + + // =============== Callbacks and Jacobian =============== + setObjEvalCallback(new CallbackEvalFC(this, completeEquationsToSolve, nonlinearConstraintIndexes)); + + List jacCstDense = new ArrayList<>(); + List jacVarDense = new ArrayList<>(); + List jacCstSparse = new ArrayList<>(); + List jacVarSparse = new ArrayList<>(); + + setJacobianMatrix( + network, jacobianMatrix, completeEquationsToSolve, nonlinearConstraintIndexes, + jacCstDense, jacVarDense, jacCstSparse, jacVarSparse + ); + + AbstractMap.SimpleEntry, List> hessNnz = getHessNnzRowsAndCols(nonlinearConstraintIndexes, completeEquationsToSolve); + setHessNnzPattern(hessNnz.getKey(), hessNnz.getValue()); + long end = System.nanoTime(); + time += end - start; + } + + /** + * Adds quadratic and linear terms related to slack variables to the objective function. + */ + private void addSlackObjectiveTerms( + int numEquations, + int slackStartIdx, + double weight, + double lambda, + double mu, + List quadRows, + List quadCols, + List quadCoefs, + List linIndexes, + List linCoefs) { + + long start = System.nanoTime(); + for (int i = 0; i < numEquations; i++) { + int idxSm = slackStartIdx + 2 * i; + int idxSp = slackStartIdx + 2 * i + 1; + + // Quadratic terms: weight * mu * (sp^2 + sm^2 - 2 * sp * sm) + quadRows.add(idxSp); + quadCols.add(idxSp); + quadCoefs.add(mu * weight); + + quadRows.add(idxSm); + quadCols.add(idxSm); + quadCoefs.add(mu * weight); + + quadRows.add(idxSp); + quadCols.add(idxSm); + quadCoefs.add(-2 * mu * weight); + + // Linear terms: weight * lambda * (sp + sm) + linIndexes.add(idxSp); + linCoefs.add(lambda * weight); + + linIndexes.add(idxSm); + linCoefs.add(lambda * weight); + } + long end = System.nanoTime(); + time += end - start; + } + + /** + * 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. + * @param equationsToSolve Base list of all equations to solve. + * @param solverUtils Utilities to extract linear constraint components. + * @param nonLinearConstraintIds Output list of non-linear constraint indices. + */ + private void addActivatedConstraints( + int equationId, + List> equationsToSolve, + NonLinearExternalSolverUtils solverUtils, + List nonLinearConstraintIds, + List> completeEquationsToSolve, + TargetVector targetVector, + List wholeTargetVector, + List listBusesWithQEqToAdd) { + + long start = System.nanoTime(); + + Equation equation = equationsToSolve.get(equationId); + AcEquationType equationType = equation.getType(); + List> terms = equation.getTerms(); + boolean addComplConstraintsVariable = listBusesWithQEqToAdd.contains(equation.getElementNum()); + + if (equationType == BUS_TARGET_V) { + 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. + + // ---- V_inf Equation ---- + if (addComplConstraintsVariable) { + int compVarBaseIndex = compVarStartIndex + 5 * vSuppEquationLocalIds.get(equationId); + varVInfIndices.add(compVarBaseIndex); // V_inf + varVInfIndices.add(compVarBaseIndex + 4); // V_aux + coefficientsVInf.add(1.0); + coefficientsVInf.add(-1.0); + } + // Add slack variables if applicable + int slackBase = getSlackIndexBase(equationType, equationId); + if (slackBase >= 0) { + varVInfIndices.add(slackBase); // Sm + varVInfIndices.add(slackBase + 1); // Sp + coefficientsVInf.add(-1.0); + coefficientsVInf.add(1.0); + } + + 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); + + // ---- V_sup Equation ---- + if (addComplConstraintsVariable) { + int compVarBaseIndex = compVarStartIndex + 5 * vSuppEquationLocalIds.get(equationId); + List varVSupIndices = new ArrayList<>(linearConstraint.listIdVar()); + List coefficientsVSup = new ArrayList<>(linearConstraint.listCoef()); + + // Add complementarity constraints' variables + 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 + if (slackBase >= 0) { + varVSupIndices.add(slackBase); // Sm + varVSupIndices.add(slackBase + 1); // Sp + coefficientsVSup.add(-1.0); + coefficientsVSup.add(1.0); + } + + for (int i = 0; i < varVSupIndices.size(); i++) { + this.addConstraintLinearPart(equationsToSolve.size() + vSuppEquationLocalIds + .get(equationId), varVSupIndices.get(i), coefficientsVSup.get(i)); + } + } + } catch (UnsupportedOperationException e) { + throw new PowsyblException("Failed to process linear constraint for equation #" + equationId, e); + } + + } else { + nonLinearConstraintIds.add(equationId); + nonLinearConstraintIds.add(equationsToSolve.size() + vSuppEquationLocalIds.get(equationId)); + } + + // Add the duplicated equation (ie V_sup eq) to the list of eq to solve and its target + if (addComplConstraintsVariable) { + Equation vEq = equationsToSolve.get(equationId); + completeEquationsToSolve.add(vEq); + wholeTargetVector.add(Arrays.stream(targetVector.getArray()).boxed().toList().get(equationId)); + } + + } else { + if (NonLinearExternalSolverUtils.isLinear(equationType, terms)) { + try { + // Extract linear constraint components + var linearConstraint = solverUtils.getLinearConstraint(equationType, terms); + List varIndices = new ArrayList<>(linearConstraint.listIdVar()); + List coefficients = new ArrayList<>(linearConstraint.listCoef()); + + // 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); + } + + for (int i = 0; i < varIndices.size(); i++) { + this.addConstraintLinearPart(equationId, varIndices.get(i), coefficients.get(i)); + } + + LOGGER.trace("Added linear constraint #{} of type {}", equationId, equationType); + } catch (UnsupportedOperationException e) { + throw new PowsyblException("Failed to process linear constraint for equation #" + equationId, e); + } + } else { + nonLinearConstraintIds.add(equationId); + } + } + long end = System.nanoTime(); + time += end - start; + } + + /** + * 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; + }; + } + + private int getcompVarBaseIndex(int equationId) { + return compVarStartIndex + 5 * vSuppEquationLocalIds.get(equationId); + } + + /** + * Configures the Jacobian matrix for the Knitro problem, using either a dense or sparse representation. + * + * @param lfNetwork The PowSyBl network. + * @param jacobianMatrix The PowSyBl Jacobian matrix. + * @param sortedEquationsToSolve The list of equations to solve. + * @param listNonLinearConsts The list of non-linear constraint ids. + * @param listNonZerosCtsDense Dense non-zero constraints. + * @param listNonZerosVarsDense Dense non-zero variables. + * @param listNonZerosCtsSparse Sparse non-zero constraints. + * @param listNonZerosVarsSparse Sparse non-zero variables. + * @throws KNException If an error occurs in Knitro operations. + */ + private void setJacobianMatrix(LfNetwork lfNetwork, JacobianMatrix jacobianMatrix, + List> sortedEquationsToSolve, List listNonLinearConsts, + List listNonZerosCtsDense, List listNonZerosVarsDense, + List listNonZerosCtsSparse, List listNonZerosVarsSparse) throws KNException { + + long start = System.nanoTime(); + + int numVar = equationSystem.getIndex().getSortedVariablesToFind().size(); + if (knitroParameters.getGradientComputationMode() == 1) { // User routine to compute the Jacobian + if (knitroParameters.getGradientUserRoutine() == 1) { + // Dense method: all non-linear constraints are considered as a function of all variables. + buildDenseJacobianMatrix(numVar, listNonLinearConsts, listNonZerosCtsDense, listNonZerosVarsDense); + this.setJacNnzPattern(listNonZerosCtsDense, listNonZerosVarsDense); + } else if (knitroParameters.getGradientUserRoutine() == 2) { + // Sparse method: compute Jacobian only for variables the constraints depend on. + buildSparseJacobianMatrix(sortedEquationsToSolve, listNonLinearConsts, listNonZerosCtsSparse, listNonZerosVarsSparse); + this.setJacNnzPattern(listNonZerosCtsSparse, listNonZerosVarsSparse); + } + // Set the callback for gradient evaluations if the user directly passes the Jacobian to the solver. + this.setGradEvalCallback(new CallbackEvalG( + jacobianMatrix, + listNonZerosCtsDense, + listNonZerosVarsDense, + listNonZerosCtsSparse, + listNonZerosVarsSparse, + lfNetwork, + equationSystem, + knitroParameters + )); + } + long end = System.nanoTime(); + time += end - start; + } + + /** + * Callback used by Knitro to evaluate the non-linear parts of the objective and constraint functions. + */ + private static final class CallbackEvalFC extends KNEvalFCCallback { + + private final List> sortedEquationsToSolve; + private final List nonLinearConstraintIds; + private final ResilientReacLimKnitroProblem problemInstance; + + private CallbackEvalFC(ResilientReacLimKnitroProblem problemInstance, List> sortedEquationsToSolve, List nonLinearConstraintIds) { + this.problemInstance = problemInstance; + this.sortedEquationsToSolve = sortedEquationsToSolve; + this.nonLinearConstraintIds = nonLinearConstraintIds; + } + + /** + * Knitro callback function that evaluates the non-linear constraints at the current point. + * + * @param x Current point (primal variables). + * @param obj Output objective value (unused here). + * @param c Output constraint values. + */ + @Override + public void evaluateFC(final List x, final List obj, final List c) { + LOGGER.trace("============ Knitro evaluating callback function ============"); + + StateVector currentState = new StateVector(toArray(x)); + LOGGER.trace("Current state vector: {}", currentState.get()); + LOGGER.trace("Evaluating {} non-linear constraints", nonLinearConstraintIds.size()); + + int callbackConstraintIndex = 0; + int nmbreEqUnactivated = sortedEquationsToSolve.stream().filter( + e -> !e.isActive()).toList().size(); + + for (int equationId : nonLinearConstraintIds) { + Equation equation = sortedEquationsToSolve.get(equationId); + AcEquationType type = equation.getType(); + + // Ensure the constraint is non-linear + if (NonLinearExternalSolverUtils.isLinear(type, equation.getTerms())) { + throw new IllegalArgumentException( + "Equation of type " + type + " is linear but passed to the non-linear callback"); + } + + // Evaluate equation using the current state + double constraintValue = 0.0; + for (EquationTerm term : equation.getTerms()) { + term.setStateVector(currentState); + if (term.isActive()) { + constraintValue += term.eval(); + } + } + + if (equation.isActive()) { // add slack variables + int slackIndexBase = problemInstance.getSlackIndexBase(type, equationId); + if (slackIndexBase >= 0) { + double sm = x.get(slackIndexBase); // negative slack + double sp = x.get(slackIndexBase + 1); // positive slack + constraintValue += sp - sm; // add s + // slack contribution + } + } else { // add blow / bup depending on the constraint + int elemNum = equation.getElementNum(); + Equation equationV = sortedEquationsToSolve.stream().filter( + e -> e.getElementNum() == elemNum).filter( + e -> e.getType() == BUS_TARGET_V).toList().get(0); + int equationVId = sortedEquationsToSolve.indexOf(equationV); + int compVarBaseIndex = problemInstance.getcompVarBaseIndex(equationVId); + if (equationId - sortedEquationsToSolve.size() + nmbreEqUnactivated / 2 < 0) { // Q_low Constraint + double bLow = x.get(compVarBaseIndex + 2); + constraintValue -= bLow; + } else { // Q_up Constraint + double bUp = x.get(compVarBaseIndex + 3); + constraintValue += bUp; + } + } + try { + c.set(callbackConstraintIndex, constraintValue); + LOGGER.trace("Added non-linear constraint #{} (type: {}) = {}", equationId, type, constraintValue); + } catch (Exception e) { + throw new PowsyblException("Error while adding non-linear constraint #" + equationId, e); + } + callbackConstraintIndex++; + } + } + } + + /** + * Callback used by Knitro to evaluate the gradient (Jacobian matrix) of the constraints. + * Only constraints (no objective) are handled here. + */ + private static final class CallbackEvalG extends KNEvalGACallback { + + private final JacobianMatrix jacobianMatrix; + private final List denseConstraintIndices; + private final List denseVariableIndices; + private final List sparseConstraintIndices; + private final List sparseVariableIndices; + + private final LfNetwork network; + private final EquationSystem equationSystem; + private final KnitroSolverParameters knitroParameters; + + private final int numLFVar; + private final int numVEq; + private final int numPQEq; + + private final LinkedHashMap indRowVariable; + + private CallbackEvalG( + JacobianMatrix jacobianMatrix, + List denseConstraintIndices, + List denseVariableIndices, + List sparseConstraintIndices, + List sparseVariableIndices, + LfNetwork network, + EquationSystem equationSystem, + KnitroSolverParameters knitroParameters) { + + this.jacobianMatrix = jacobianMatrix; + this.denseConstraintIndices = denseConstraintIndices; + this.denseVariableIndices = denseVariableIndices; + this.sparseConstraintIndices = sparseConstraintIndices; + this.sparseVariableIndices = sparseVariableIndices; + this.network = network; + this.equationSystem = equationSystem; + this.knitroParameters = knitroParameters; + + 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(); + + this.indRowVariable = new LinkedHashMap(); + List> listVar = equationSystem.getVariableSet().getVariables().stream().toList(); + for (Variable variable : listVar) { + indRowVariable.put(variable.getRow(), variable); + } + } + + @Override + public void evaluateGA(final List x, final List objGrad, final List jac) { + // Update internal state and Jacobian + equationSystem.getStateVector().set(toArray(x)); + AcSolverUtil.updateNetwork(network, equationSystem); + jacobianMatrix.forceUpdate(); + + // Get sparse matrix representation + SparseMatrix sparseMatrix = jacobianMatrix.getMatrix().toSparse(); + int[] columnStart = sparseMatrix.getColumnStart(); + double[] values = sparseMatrix.getValues(); + + // Determine which list to use based on Knitro settings + List constraintIndices; + List variableIndices; + + int routineType = knitroParameters.getGradientUserRoutine(); + if (routineType == 1) { + constraintIndices = denseConstraintIndices; + variableIndices = denseVariableIndices; + } else if (routineType == 2) { + constraintIndices = sparseConstraintIndices; + variableIndices = sparseVariableIndices; + } else { + throw new IllegalArgumentException("Unsupported gradientUserRoutine value: " + routineType); + } + + // Fill Jacobian values + boolean firstIteration = true; + int iRowIndices = 0; + int currentConstraint = -1; + int numVariables = equationSystem.getIndex().getSortedVariablesToFind().size(); + for (int index = 0; index < constraintIndices.size(); index++) { + try { + if (firstIteration) { + currentConstraint = constraintIndices.get(index); + } + + int ct = constraintIndices.get(index); + int var = variableIndices.get(index); + double value; + + if (ct < Arrays.stream(columnStart).count()) { + + // Find matching (var, ct) entry in sparse column + int colStart = columnStart[ct]; + + if (!firstIteration) { + if (currentConstraint != ct) { + iRowIndices = 0; + currentConstraint = ct; + } + } + + if (var >= numVariables) { + if (((var - numVariables) & 1) == 0) { + // set Jacobian entry to -1.0 if slack variable is Sm + value = -1.0; + } else { + // set Jacobian entry to +1.0 if slack variable is Sp + value = 1.0; + } + } else { + value = values[colStart + iRowIndices++]; + } + jac.set(index, value); + if (firstIteration) { + firstIteration = false; + } + } else { + // If var is a LF variable : derivate non-activated equations + value = 0.0; + Equation equation = INDEQUNACTIVEQ.get(ct); + if (var < 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(var) == v) { + //equationSystem.getVariableSet().getVariables().stream().filter(va -> va.getRow() == var ).toList().get(0) == v) { + value += term.isActive() ? term.der(v) : 0; + } + + } + } + } + // Check if var is a b_low or b_up var + if (var >= numLFVar + 2 * (numPQEq + numVEq)) { + int rest = (var - 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; + } + } + jac.set(index, value); + } + } catch (Exception e) { + int varId = routineType == 1 ? denseVariableIndices.get(index) : sparseVariableIndices.get(index); + int ctId = routineType == 1 ? denseConstraintIndices.get(index) : sparseConstraintIndices.get(index); + LOGGER.error("Error while filling Jacobian term at var {} and constraint {}", varId, ctId, e); + } + } + } + } + } +} diff --git a/src/main/java/com/powsybl/openloadflow/knitro/solver/KnitroWritter.java b/src/main/java/com/powsybl/openloadflow/knitro/solver/KnitroWritter.java new file mode 100644 index 0000000..53ebc78 --- /dev/null +++ b/src/main/java/com/powsybl/openloadflow/knitro/solver/KnitroWritter.java @@ -0,0 +1,26 @@ +package com.powsybl.openloadflow.knitro.solver; + +import java.io.BufferedWriter; +import java.io.FileWriter; +import java.io.IOException; + +public class KnitroWritter { + private final String logFile; + + public KnitroWritter(String logFile) { + this.logFile = logFile; + } + + public void write(String message, boolean append) { + try (BufferedWriter writer = new BufferedWriter(new FileWriter(logFile, append))) { + writer.write(message); + writer.newLine(); + } catch (IOException e) { + System.err.println("Erreur lors de l'écriture des logs : " + e.getMessage()); + } + } + + public String getLogFile() { + return logFile; + } +} diff --git a/src/main/java/com/powsybl/openloadflow/knitro/solver/ReacLimitsTestsUtils.java b/src/main/java/com/powsybl/openloadflow/knitro/solver/ReacLimitsTestsUtils.java new file mode 100644 index 0000000..b7c6a27 --- /dev/null +++ b/src/main/java/com/powsybl/openloadflow/knitro/solver/ReacLimitsTestsUtils.java @@ -0,0 +1,233 @@ +/** + * 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.iidm.network.*; +import com.powsybl.loadflow.LoadFlow; +import com.powsybl.loadflow.LoadFlowParameters; +import com.powsybl.loadflow.LoadFlowResult; +import com.powsybl.openloadflow.OpenLoadFlowParameters; +import com.powsybl.openloadflow.network.PlausibleValues; +import org.slf4j.Logger; +import org.slf4j.LoggerFactory; + +import java.io.*; +import java.util.*; + + +import static org.ejml.UtilEjml.assertTrue; + +/** + * @author Pierre Arvy {@literal } + * @author Yoann Anezin {@literal } + */ +public final class ReacLimitsTestsUtils { + private static final double DEFAULT_TOLERANCE = 1e-4; + private static final Logger LOGGER = LoggerFactory.getLogger(ReacLimitsTestsUtils.class); + + private ReacLimitsTestsUtils() { + throw new UnsupportedOperationException(); + } + + public static ArrayList countAndSwitch(Network network, HashMap listMinQ, HashMap listMaxQ, + HashMap slacksP, HashMap slacksQ, HashMap slacksV) throws Exception { + int nmbSwitchQmin = 0; + int nmbSwitchQmax = 0; + int previousNmbBusPV = 0; + ArrayList switches = new ArrayList<>(); + HashMap visitedBuses = new HashMap<>(); + for (Generator g : network.getGenerators()) { + if (g.getTerminal().getBusView().getBus() == null || !g.isVoltageRegulatorOn()) { + continue; + } + String idbus = g.getTerminal().getBusView().getBus().getId(); + Double slackP = slacksP.get(idbus); + Double slackQ = slacksQ.get(idbus); + Double slackV = slacksV.get(idbus); + if (slackP == null) { + slackP = 0.0; + } + if (slackQ == null) { + slackQ = 0.0; + } + if (slackV == null) { + slackV = 0.0; + } + Terminal t = g.getTerminal(); + double v = t.getBusView().getBus().getV(); + double qMing = g.getReactiveLimits().getMinQ(g.getTargetP()); + double qMaxg = g.getReactiveLimits().getMaxQ(g.getTargetP()); + if (Math.abs(qMing - qMaxg) < PlausibleValues.MIN_REACTIVE_RANGE) { + continue; + } + if (Math.abs(qMing) > PlausibleValues.MAX_REACTIVE_RANGE || Math.abs(qMaxg) > PlausibleValues.MAX_REACTIVE_RANGE) { + continue; + } + + if (visitedBuses.containsKey(idbus)) { + switch (visitedBuses.get(idbus)) { + case "Switch Qmin": + assertTrue(v + slackV > g.getTargetV(), "Another generator did a Qmin switch," + + " expected the same thing to happened. Current generator : " + g.getId() + " on bus " + idbus); + assertTrue(-t.getQ() + DEFAULT_TOLERANCE > qMing && -t.getQ() - DEFAULT_TOLERANCE < qMing, + "Another generator did a Qmin switch, expected the same thing to happened. " + + "Current generator : " + g.getId() + " on bus " + idbus); + break; + case "Switch Qmax": + assertTrue(v + slackV < g.getTargetV(), "Another generator did a Qmax switch," + + " expected the same thing to happened. Current generator : " + g.getId() + " on bus " + idbus); + assertTrue(-t.getQ() + DEFAULT_TOLERANCE > qMaxg && -t.getQ() - DEFAULT_TOLERANCE < qMaxg, + "Another generator did a Qmax switch, expected the same thing to happened. " + + "Current generator : " + g.getId() + " on bus " + idbus); + break; + case "No Switch": + assertTrue(v + slackV + DEFAULT_TOLERANCE > g.getTargetV() && v + slackV - DEFAULT_TOLERANCE < g.getTargetV()); + } + continue; + } + previousNmbBusPV++; + if (!(v + slackV + DEFAULT_TOLERANCE > g.getTargetV() && v + slackV - DEFAULT_TOLERANCE < g.getTargetV())) { + + if ((-t.getQ() + 2*DEFAULT_TOLERANCE > qMing && + -t.getQ() - 2*DEFAULT_TOLERANCE < qMing) && qMing != qMaxg) { + nmbSwitchQmin++; + if (!(v + slackV > g.getTargetV())) { + LOGGER.warn("V ( " + v + " ) below its target ( " + g.getTargetV() + " ) on a Qmin switch of bus " + + t.getBusView().getBus().getId() + ". Current generator checked : " + g.getId()); + } + g.setTargetQ(listMinQ.get(g.getId())); + visitedBuses.put(idbus, "Switch Qmin"); + } else if ((-t.getQ() + DEFAULT_TOLERANCE > qMaxg && + -t.getQ() - DEFAULT_TOLERANCE < qMaxg) && qMing != qMaxg) { + nmbSwitchQmax++; + if (!(v + slackV < g.getTargetV())) { + LOGGER.warn("V ( " + v + " ) above its target ( " + g.getTargetV() + " ) on a Qmax switch of bus " + + t.getBusView().getBus().getId() + ". Current generator checked : " + g.getId()); + } + g.setTargetQ(listMaxQ.get(g.getId())); + visitedBuses.put(idbus, "Switch Qmax"); + } else if ((-t.getQ() + DEFAULT_TOLERANCE > qMaxg && + -t.getQ() - DEFAULT_TOLERANCE < qMaxg) && qMaxg == qMing) { + if (v + slackV > g.getTargetV()) { + nmbSwitchQmin++; + g.setTargetQ(listMinQ.get(g.getId())); + visitedBuses.put(idbus, "Switch Qmin"); + } else { + nmbSwitchQmax++; + g.setTargetQ(listMaxQ.get(g.getId())); + visitedBuses.put(idbus, "Switch Qmax"); + } + } else { + assertTrue((-t.getQ() + DEFAULT_TOLERANCE > qMaxg && -t.getQ() - DEFAULT_TOLERANCE < qMaxg) || + (-t.getQ() + DEFAULT_TOLERANCE > qMing && -t.getQ() - DEFAULT_TOLERANCE < qMing), + "Value of Q ( " + -t.getQ() + " ) not matching Qmin ( " + qMing + " ) nor Qmax ( " + + qMaxg + " ) on the switch of bus " + t.getBusView().getBus().getId() + + ". Current generator checked : " + g.getId()); + } + g.setVoltageRegulatorOn(false); + } else { + visitedBuses.put(idbus, "No Switch"); + } + } + switches.add(nmbSwitchQmin); + switches.add(nmbSwitchQmax); + switches.add(previousNmbBusPV); + return switches; + } + + public static void verifNewtonRaphson(Network network, LoadFlowParameters parameters, LoadFlow.Runner loadFlowRunner, int nbreIter) { + OpenLoadFlowParameters.get(parameters).setVoltageInitModeOverride(OpenLoadFlowParameters.VoltageInitModeOverride.NONE); + parameters.setVoltageInitMode(LoadFlowParameters.VoltageInitMode.PREVIOUS_VALUES); + OpenLoadFlowParameters.get(parameters).setMaxNewtonRaphsonIterations(nbreIter) + .setReportedFeatures(Collections.singleton(OpenLoadFlowParameters.ReportedFeatures.NEWTON_RAPHSON_LOAD_FLOW)); + OpenLoadFlowParameters.get(parameters).setAcSolverType("NEWTON_RAPHSON"); + LoadFlowResult result = loadFlowRunner.run(network, parameters); + assertTrue(result.isFullyConverged()); + } + + public static void checkSwitches(Network network, HashMap listMinQ, HashMap listMaxQ) { + HashMap slacksP = new HashMap(); + HashMap slacksQ = new HashMap(); + HashMap slacksV = new HashMap(); + ArrayList slacksfiles = new ArrayList(); + slacksfiles.add("P"); + slacksfiles.add("Q"); + slacksfiles.add("V"); + for (String type : slacksfiles) { + try (BufferedReader br = new BufferedReader(new FileReader("D:\\Documents\\Slacks\\Slacks" + type + ".txt"))) { + String ligne; + boolean isId = true; // True = on attend un identifiant, False = on attend une valeur + String id = ""; + Double value; + + while ((ligne = br.readLine()) != null) { + if (ligne.trim().isEmpty()) { + continue; // ignorer les lignes vides éventuelles + } + + if (isId) { + id = ligne.trim(); + } else { + // Convertir la valeur en double + try { + value = Double.parseDouble(ligne.trim().replace(',', '.')); + switch (type) { + case "P": + slacksP.put(id, value); + break; + case "Q": + slacksQ.put(id, value); + break; + case "V": + slacksV.put(id, value); + break; + } + } catch (NumberFormatException e) { + System.err.println("Valeur non numérique : " + ligne); + } + } + isId = !isId; // alterner ID/valeur + } + + } catch (IOException e) { + e.printStackTrace(); + } + } + + try { + ArrayList switches = countAndSwitch(network, listMinQ, listMaxQ, slacksP, slacksQ, slacksV); + assertTrue(switches.get(2) > switches.get(1) + switches.get(0), + "No control on any voltage magnitude : all buses switched"); + System.out.println(switches.get(0) + " switches to PQ with Q = Qlow and " + switches.get(1) + " with Q = Qup"); + } catch (Exception e) { + throw new RuntimeException(e); + } + + for (String type : slacksfiles) { + try (FileWriter fw = new FileWriter("D:\\Documents\\Slacks\\Slacks" + type + ".txt", false)) { + fw.write(""); + } catch (IOException e) { + e.printStackTrace(); + } + } + } + + /** + * Creates an active power perturbation of a given network. + * + * @param network The network to perturb. + * @param alpha The active load mismatch to apply. + */ + public static void applyActivePowerPerturbation(Network network, double alpha) { + for (Load load : network.getLoads()) { + load.setP0(alpha * load.getP0()); + } + } +} diff --git a/src/main/java/com/powsybl/openloadflow/knitro/solver/ResilientKnitroSolver.java b/src/main/java/com/powsybl/openloadflow/knitro/solver/ResilientKnitroSolver.java index 5726e73..92edff5 100644 --- a/src/main/java/com/powsybl/openloadflow/knitro/solver/ResilientKnitroSolver.java +++ b/src/main/java/com/powsybl/openloadflow/knitro/solver/ResilientKnitroSolver.java @@ -53,7 +53,7 @@ public class ResilientKnitroSolver extends AbstractAcSolver { private final double wK = 1.0; private final double wP = 1.0; private final double wQ = 1.0; - private final double wV = 1.0; + private final double wV = 10.0; // Weights of the linear and quadratic terms in the objective function private final double lambda = 3.0; @@ -115,9 +115,9 @@ public ResilientKnitroSolver( this.slackVStartIndex = slackQStartIndex + 2 * numQEquations; // Map equations to local indices - this.pEquationLocalIds = new HashMap<>(); - this.qEquationLocalIds = new HashMap<>(); - this.vEquationLocalIds = new HashMap<>(); + this.pEquationLocalIds = new LinkedHashMap<>(); + this.qEquationLocalIds = new LinkedHashMap<>(); + this.vEquationLocalIds = new LinkedHashMap<>(); int pCounter = 0; int qCounter = 0; @@ -243,12 +243,15 @@ private void setSolverParameters(KNSolver solver) throws KNException { solver.setParam(KNConstants.KN_PARAM_FEASTOL, knitroParameters.getConvEps()); solver.setParam(KNConstants.KN_PARAM_MAXIT, knitroParameters.getMaxIterations()); solver.setParam(KNConstants.KN_PARAM_HESSOPT, knitroParameters.getHessianComputationMode()); - solver.setParam(KNConstants.KN_PARAM_SOLTYPE, KNConstants.KN_SOLTYPE_BESTFEAS); +// solver.setParam(KNConstants.KN_PARAM_SOLTYPE, KNConstants.KN_SOLTYPE_BESTFEAS); solver.setParam(KNConstants.KN_PARAM_OUTMODE, KNConstants.KN_OUTMODE_BOTH); solver.setParam(KNConstants.KN_PARAM_OPTTOL, 1.0e-3); solver.setParam(KNConstants.KN_PARAM_OPTTOLABS, 1.0e-1); solver.setParam(KNConstants.KN_PARAM_OUTLEV, 3); - solver.setParam(KNConstants.KN_PARAM_NUMTHREADS, 8); + solver.setParam(KNConstants.KN_PARAM_NUMTHREADS, 1); +// solver.setParam(KNConstants.KN_PARAM_ALG, 1); +// solver.setParam(KNConstants.KN_PARAM_DERIVCHECK, 1); +// solver.setParam(KNConstants.KN_PARAM_FINDIFF_NUMTHREADS, 1); LOGGER.info("Knitro parameters set: GRADOPT={}, HESSOPT={}, FEASTOL={}, MAXIT={}", knitroParameters.getGradientComputationMode(), @@ -274,7 +277,7 @@ public AcSolverResult run(VoltageInitializer voltageInitializer, ReportNode repo setSolverParameters(solver); solver.solve(); - KNSolution solution = solver.getBestFeasibleIterate(); + KNSolution solution = solver.getSolution(); List x = solution.getX(); solverStatus = KnitroStatus.fromStatusCode(solution.getStatus()).toAcSolverStatus(); @@ -302,7 +305,6 @@ public AcSolverResult run(VoltageInitializer voltageInitializer, ReportNode repo LOGGER.info("Penalty Q = {}", penaltyQ); LOGGER.info("Penalty V = {}", penaltyV); LOGGER.info("Total penalty = {}", totalPenalty); - // ========== Network Update ========== // Update the network values if the solver converged or if the network should always be updated if (solverStatus == AcSolverStatus.CONVERGED || knitroParameters.isAlwaysUpdateNetwork()) { diff --git a/src/test/java/com/powsybl/openloadflow/knitro/solver/AcloadFlowReactiveLimitsTest.java b/src/test/java/com/powsybl/openloadflow/knitro/solver/AcloadFlowReactiveLimitsTest.java index 221b1f3..062007c 100644 --- a/src/test/java/com/powsybl/openloadflow/knitro/solver/AcloadFlowReactiveLimitsTest.java +++ b/src/test/java/com/powsybl/openloadflow/knitro/solver/AcloadFlowReactiveLimitsTest.java @@ -107,6 +107,7 @@ void setUp() { .setDistributedSlack(false); KnitroLoadFlowParameters knitroLoadFlowParameters = new KnitroLoadFlowParameters(); // set gradient computation mode knitroLoadFlowParameters.setGradientComputationMode(2); + knitroLoadFlowParameters.setKnitroSolverType(KnitroSolverParameters.KnitroSolverType.REACTIVLIMITS); parameters.addExtension(KnitroLoadFlowParameters.class, knitroLoadFlowParameters); OpenLoadFlowParameters.create(parameters).setAcSolverType(KnitroSolverFactory.NAME); } @@ -125,15 +126,27 @@ void diagramTest() { @Test void test() { - parameters.setUseReactiveLimits(false); + + /*parameters.setUseReactiveLimits(false); + LoadFlowResult result = loadFlowRunner.run(network, parameters); + assertTrue(result.isFullyConverged()); + assertReactivePowerEquals(-109.229, gen.getTerminal()); + assertReactivePowerEquals(-152.266, gen2.getTerminal()); + assertReactivePowerEquals(-199.999, nhv2Nload.getTerminal2());*/ + + /*parameters.setUseReactiveLimits(true); + parameters.getExtension(OpenLoadFlowParameters.class) + .setAlwaysUpdateNetwork(true) + .setAcSolverType(NewtonRaphsonFactory.NAME); LoadFlowResult result = loadFlowRunner.run(network, parameters); assertTrue(result.isFullyConverged()); - assertReactivePowerEquals(-109.228, gen.getTerminal()); - assertReactivePowerEquals(-152.265, gen2.getTerminal()); - assertReactivePowerEquals(-199.998, nhv2Nload.getTerminal2()); + parameters.setVoltageInitMode(LoadFlowParameters.VoltageInitMode.PREVIOUS_VALUES); + parameters.getExtension(OpenLoadFlowParameters.class) + .setAlwaysUpdateNetwork(true) + .setAcSolverType(KnitroSolverFactory.NAME);*/ parameters.setUseReactiveLimits(true); - result = loadFlowRunner.run(network, parameters); + LoadFlowResult result = loadFlowRunner.run(network, parameters); assertTrue(result.isFullyConverged()); assertReactivePowerEquals(-164.315, gen.getTerminal()); assertReactivePowerEquals(-100, gen2.getTerminal()); // GEN is correctly limited to 100 MVar @@ -163,4 +176,9 @@ void testWithMixedGenLoad() { assertReactivePowerEquals(-120, gen2.getTerminal()); assertReactivePowerEquals(100, ngen2Nhv1.getTerminal1()); } + + private boolean checkNotOnlyPQ(LfNetwork lfNetwork, LoadFlowResult result) { + + return false; + } } diff --git a/src/test/java/com/powsybl/openloadflow/knitro/solver/ReacLimPertubationTest.java b/src/test/java/com/powsybl/openloadflow/knitro/solver/ReacLimPertubationTest.java new file mode 100644 index 0000000..a41d1ed --- /dev/null +++ b/src/test/java/com/powsybl/openloadflow/knitro/solver/ReacLimPertubationTest.java @@ -0,0 +1,625 @@ +/** + * 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.loadflow.LoadFlow; +import com.powsybl.loadflow.LoadFlowParameters; +import com.powsybl.loadflow.LoadFlowResult; +import com.powsybl.math.matrix.SparseMatrixFactory; +import com.powsybl.openloadflow.OpenLoadFlowParameters; +import com.powsybl.openloadflow.OpenLoadFlowProvider; +import com.powsybl.openloadflow.network.SlackBusSelectionMode; +import org.junit.jupiter.api.BeforeEach; +import org.junit.jupiter.api.Test; + +import java.io.IOException; +import java.time.LocalDateTime; + +import java.util.*; + +import static org.junit.jupiter.api.Assertions.*; + +/** + * @author Pierre Arvy {@literal } + * @author Yoann Anezin {@literal } + */ +public class ReacLimPertubationTest { + private LoadFlow.Runner loadFlowRunner; + private LoadFlowParameters parameters; + private static String logFile = "D:\\Documents\\Logs_Tests\\Logs.txt"; + + public ReacLimPertubationTest() throws IOException { + } + + private int fixReacLim(Network network, HashMap listMinQ, HashMap listMaxQ) { + int numbreLimReacAdded = 0; + for (var g : network.getGenerators()) { + if (g.getReactiveLimits().getMinQ(g.getTargetP()) > -1.7976931348623157E308) { + listMinQ.put(g.getId(), g.getReactiveLimits().getMinQ(g.getTargetP())); + listMaxQ.put(g.getId(), g.getReactiveLimits().getMaxQ(g.getTargetP())); + } else { + g.newMinMaxReactiveLimits() + .setMinQ(-2000) + .setMaxQ(2000) + .add(); + listMinQ.put(g.getId(), -2000.0); + listMaxQ.put(g.getId(), 2000.0); + numbreLimReacAdded++; + } + } + return numbreLimReacAdded; + } + + void logsWriting(KnitroLoadFlowParameters knitroLoadFlowParameters, KnitroWritter knitroWritter) { + knitroWritter.write("Solver used : " + OpenLoadFlowParameters.get(parameters).getAcSolverType(), true); + knitroWritter.write("Init Mode : " + parameters.getVoltageInitMode().toString(), true); + if (OpenLoadFlowParameters.get(parameters).getAcSolverType().equals("KNITRO")) { + knitroWritter.write("Version de Knitro : " + knitroLoadFlowParameters.getKnitroSolverType().name(), true); + } + knitroWritter.write("Override Init Mode : " + OpenLoadFlowParameters.get(parameters).getVoltageInitModeOverride().name(), true); + knitroWritter.write("Max Iterations : " + knitroLoadFlowParameters.getMaxIterations(), true); + knitroWritter.write("Gradient Computation Mode : " + knitroLoadFlowParameters.getGradientComputationMode(), true); + } + + /** + *Start all the test process and writes logs by the same time + * @param logFile file where logs are written + * @param network network of work + * @param perturbProcess Indicates the perturbation to apply. Current possible choices : ActivePowerGlobal, + * ActivePowerLocal, ReactivePower, None + * @param perturbValue value applied in the pertubation chosen (wont be used in the "None" case) + */ + void testprocess(String logFile, Network network, String perturbProcess, double perturbValue) { + long start = System.nanoTime(); + + KnitroWritter knitroWritter = new KnitroWritter(logFile); + KnitroLoadFlowParameters knitroLoadFlowParameters = parameters.getExtension(KnitroLoadFlowParameters.class); + knitroLoadFlowParameters.setKnitroWritter(knitroWritter); + parameters.addExtension(KnitroLoadFlowParameters.class, knitroLoadFlowParameters); + + HashMap listMinQ = new HashMap<>(); + HashMap listMaxQ = new HashMap<>(); + parameters.setUseReactiveLimits(true); + int numbreLimReacAdded = fixReacLim(network, listMinQ, listMaxQ); + + knitroWritter.write("[" + LocalDateTime.now() + "]", false); + knitroWritter.write(numbreLimReacAdded + " bus pour lesquels les limites réactif ont été ajoutées", true); + logsWriting(knitroLoadFlowParameters, knitroWritter); + switch (perturbProcess) { + case "ActivePowerGlobal": + ReacLimitsTestsUtils.applyActivePowerPerturbation(network, perturbValue); + knitroWritter.write("Perturbed by global loads' increasement (All multiplied by " + + perturbValue * 100 + "% of total load)", true); + break; + case "ActivePowerLocal": + PerturbationFactory.applyActivePowerPerturbation(network, + PerturbationFactory.getActivePowerPerturbation(network), perturbValue); + knitroWritter.write("Perturbed by uniq big load (" + perturbValue * 100 + "% of total load)", true); + break; + case "ReactivePower": + PerturbationFactory.applyReactivePowerPerturbation(network, + PerturbationFactory.getReactivePowerPerturbation(network), perturbValue); + knitroWritter.write("Perturbed by power injection by the shunt (Target Q = " + perturbValue + ")", true); + break; + case "None": + knitroWritter.write("No Pertubations", true); + break; + default: + knitroWritter.write("No Pertubations, you have miswritten the pertubation's instruction", true); + break; + } + + // solve and check + LoadFlowResult result = loadFlowRunner.run(network, parameters); + assertTrue(result.isFullyConverged(), "Not Fully Converged"); + ReacLimitsTestsUtils.checkSwitches(network, listMinQ, listMaxQ); + long end = System.nanoTime(); + knitroWritter.write("Durée du test : " + (end - start) * 1e-9 + " secondes", true); + knitroWritter.write("Nombre d'itérations : " + result.getComponentResults().get(0).getIterationCount(), true); + knitroWritter.write("Status à l'arrivée : " + result.getComponentResults().get(0).getStatus().name(), true); + } + + @BeforeEach + void setUp() { + loadFlowRunner = new LoadFlow.Runner(new OpenLoadFlowProvider(new SparseMatrixFactory())); + parameters = new LoadFlowParameters().setUseReactiveLimits(true) + .setDistributedSlack(false); + KnitroLoadFlowParameters knitroLoadFlowParameters = new KnitroLoadFlowParameters(); // set gradient computation mode + knitroLoadFlowParameters.setGradientComputationMode(1); + knitroLoadFlowParameters.setMaxIterations(2000); + knitroLoadFlowParameters.setKnitroSolverType(KnitroSolverParameters.KnitroSolverType.REACTIVLIMITS); + parameters.addExtension(KnitroLoadFlowParameters.class, knitroLoadFlowParameters); + //parameters.setVoltageInitMode(LoadFlowParameters.VoltageInitMode.DC_VALUES); + //OpenLoadFlowParameters.create(parameters).setAcSolverType("NEWTON_RAPHSON"); +// OpenLoadFlowParameters.create(parameters).setAcSolverType(KnitroSolverFactory.NAME); + OpenLoadFlowParameters.get(parameters).setVoltageInitModeOverride(OpenLoadFlowParameters.VoltageInitModeOverride.FULL_VOLTAGE); + + } + + /** + *
+     *     G1        LD2        G3
+     *     |    L12   |   L23   |
+     *     |  ------- | ------- |
+     *     B1         B2        B3
+     *
+ */ + @Test + public void createNetworkWithT2wtActivePower() { + + Network network = Network.create("yoann-n", "test"); + + Substation substation1 = network.newSubstation() + .setId("SUBSTATION1") + .setCountry(Country.FR) + .add(); + VoltageLevel vl1 = substation1.newVoltageLevel() + .setId("VL_1") + .setNominalV(132.0) + .setLowVoltageLimit(118.8) + .setHighVoltageLimit(145.2) + .setTopologyKind(TopologyKind.BUS_BREAKER) + .add(); + vl1.getBusBreakerView().newBus() + .setId("BUS_1") + .add(); + Generator g1 = vl1.newGenerator() + .setId("GEN_1") + .setBus("BUS_1") + .setMinP(0.0) + .setMaxP(140) + .setTargetP(25) + .setTargetV(135) + .setVoltageRegulatorOn(true) + .add(); + g1.newMinMaxReactiveLimits().setMinQ(-30000.0).setMaxQ(30000.0).add(); + + Substation substation = network.newSubstation() + .setId("SUBSTATION") + .setCountry(Country.FR) + .add(); + VoltageLevel vl2 = substation.newVoltageLevel() + .setId("VL_2") + .setNominalV(132.0) + .setLowVoltageLimit(118.8) + .setHighVoltageLimit(145.2) + .setTopologyKind(TopologyKind.BUS_BREAKER) + .add(); + vl2.getBusBreakerView().newBus() + .setId("BUS_2") + .add(); + vl2.newLoad() + .setId("LOAD_2") + .setBus("BUS_2") + .setP0(35) + .setQ0(20) + .add(); + + VoltageLevel vl3 = substation.newVoltageLevel() + .setId("VL_3") + .setNominalV(132.0) + .setLowVoltageLimit(118.8) + .setHighVoltageLimit(145.2) + .setTopologyKind(TopologyKind.BUS_BREAKER) + .add(); + vl3.getBusBreakerView().newBus() + .setId("BUS_3") + .add(); + Generator g3 = vl3.newGenerator() + .setId("GEN_3") + .setBus("BUS_3") + .setMinP(0.0) + .setMaxP(140) + .setTargetP(15) + .setTargetV(130) + .setVoltageRegulatorOn(true) + .add(); + g3.newMinMaxReactiveLimits().setMinQ(-3000.0).setMaxQ(3000.0).add(); + + network.newLine() + .setId("LINE_12") + .setBus1("BUS_1") + .setBus2("BUS_2") + .setR(1.05) + .setX(10.0) + .setG1(0.0000005) + .add(); + network.newLine() + .setId("LINE_23") + .setBus1("BUS_2") + .setBus2("BUS_3") + .setR(1.05) + .setX(10.0) + .setG1(0.0000005) + .add(); + + OpenLoadFlowParameters.get(parameters) +// .setVoltageInitModeOverride(OpenLoadFlowParameters.VoltageInitModeOverride.FULL_VOLTAGE) + .setSlackBusSelectionMode(SlackBusSelectionMode.NAME) + .setSlackBusId("VL_1_0"); + + logFile = "D:\\Documents\\Logs_Tests\\Logs_3bus_Active_Power_Perturbation.txt"; + testprocess(logFile, network, "ActivePowerLocal", 40.0); + ReacLimitsTestsUtils.verifNewtonRaphson(network, parameters, loadFlowRunner, 20); + } + + /** + *
+     *     G1        LD2        G3
+     *     |    L12   |   L23   |
+     *     |  ------- | ------- |
+     *     B1         B2        B3
+     *
+ */ + @Test + public void createNetworkWithT2wtReactivePower() { + HashMap listMinQ = new HashMap<>(); + HashMap listMaxQ = new HashMap<>(); + + Network network = Network.create("yoann-n", "test"); + + Substation substation1 = network.newSubstation() + .setId("SUBSTATION1") + .setCountry(Country.FR) + .add(); + VoltageLevel vl1 = substation1.newVoltageLevel() + .setId("VL_1") + .setNominalV(132.0) + .setLowVoltageLimit(118.8) + .setHighVoltageLimit(145.2) + .setTopologyKind(TopologyKind.BUS_BREAKER) + .add(); + vl1.getBusBreakerView().newBus() + .setId("BUS_1") + .add(); + Generator g1 = vl1.newGenerator() + .setId("GEN_1") + .setBus("BUS_1") + .setMinP(0.0) + .setMaxP(140) + .setTargetP(25) + .setTargetV(135) + .setVoltageRegulatorOn(true) + .add(); + g1.newMinMaxReactiveLimits().setMinQ(-3000.0).setMaxQ(3000.0).add(); + listMinQ.put(g1.getId(), -3000.0); + listMaxQ.put(g1.getId(), 3000.0); + + Substation substation = network.newSubstation() + .setId("SUBSTATION") + .setCountry(Country.FR) + .add(); + VoltageLevel vl2 = substation.newVoltageLevel() + .setId("VL_2") + .setNominalV(132.0) + .setLowVoltageLimit(118.8) + .setHighVoltageLimit(145.2) + .setTopologyKind(TopologyKind.BUS_BREAKER) + .add(); + vl2.getBusBreakerView().newBus() + .setId("BUS_2") + .add(); + vl2.newLoad() + .setId("LOAD_2") + .setBus("BUS_2") + .setP0(35) + .setQ0(20) + .add(); + + VoltageLevel vl3 = substation.newVoltageLevel() + .setId("VL_3") + .setNominalV(132.0) + .setLowVoltageLimit(118.8) + .setHighVoltageLimit(145.2) + .setTopologyKind(TopologyKind.BUS_BREAKER) + .add(); + vl3.getBusBreakerView().newBus() + .setId("BUS_3") + .add(); + Generator g3 = vl3.newGenerator() + .setId("GEN_3") + .setBus("BUS_3") + .setMinP(0.0) + .setMaxP(140) + .setTargetP(15) + .setTargetV(130) + .setVoltageRegulatorOn(true) + .add(); + g3.newMinMaxReactiveLimits().setMinQ(-3000.0).setMaxQ(3000.0).add(); + listMinQ.put(g3.getId(), -3000.0); + listMaxQ.put(g3.getId(), 3000.0); + + network.newLine() + .setId("LINE_12") + .setBus1("BUS_1") + .setBus2("BUS_2") + .setR(1.05) + .setX(10.0) + .setG1(0.0000005) + .add(); + network.newLine() + .setId("LINE_23") + .setBus1("BUS_2") + .setBus2("BUS_3") + .setR(1.05) + .setX(10.0) + .setG1(0.0000005) + .add(); + network.getVoltageLevel("VL_1").newShuntCompensator() + .setId("SC") + .setBus("BUS_1") + .setConnectableBus("BUS_1") + .setSectionCount(1) + .newLinearModel() + .setBPerSection(3.25 * Math.pow(10, -3)) + .setMaximumSectionCount(1) + .add() + .add(); + + OpenLoadFlowParameters.get(parameters) +// .setVoltageInitModeOverride(OpenLoadFlowParameters.VoltageInitModeOverride.FULL_VOLTAGE) + .setSlackBusSelectionMode(SlackBusSelectionMode.NAME) + .setSlackBusId("VL_1_0"); + + logFile = "D:\\Documents\\Logs_Tests\\Logs_3bus_Reactive_Power_Perturbation.txt"; + testprocess(logFile, network, "ReactivePower", 1E10); + ReacLimitsTestsUtils.verifNewtonRaphson(network, parameters, loadFlowRunner, 20); + } + + /** + *
+     *     G1        LD2        G3
+     *     |    L12   |   L23   |
+     *     |  ------- | ------- |
+     *     B1         B2        B3
+     *
+ */ + @Test + public void createNetworkWithT2wtVoltage() { + HashMap listMinQ = new HashMap<>(); + HashMap listMaxQ = new HashMap<>(); + Network network = Network.create("yoann-n", "test"); + + Substation substation1 = network.newSubstation() + .setId("SUBSTATION1") + .setCountry(Country.FR) + .add(); + VoltageLevel vl1 = substation1.newVoltageLevel() + .setId("VL_1") + .setNominalV(132.0) + .setLowVoltageLimit(118.8) + .setHighVoltageLimit(145.2) + .setTopologyKind(TopologyKind.BUS_BREAKER) + .add(); + vl1.getBusBreakerView().newBus() + .setId("BUS_1") + .add(); + Generator g1 = vl1.newGenerator() + .setId("GEN_1") + .setBus("BUS_1") + .setMinP(0.0) + .setMaxP(140) + .setTargetP(25) + .setTargetV(135) + .setVoltageRegulatorOn(true) + .add(); + g1.newMinMaxReactiveLimits().setMinQ(-3000.0).setMaxQ(3000.0).add(); + + Substation substation = network.newSubstation() + .setId("SUBSTATION") + .setCountry(Country.FR) + .add(); + VoltageLevel vl2 = substation.newVoltageLevel() + .setId("VL_2") + .setNominalV(132.0) + .setLowVoltageLimit(118.8) + .setHighVoltageLimit(145.2) + .setTopologyKind(TopologyKind.BUS_BREAKER) + .add(); + vl2.getBusBreakerView().newBus() + .setId("BUS_2") + .add(); + vl2.newLoad() + .setId("LOAD_2") + .setBus("BUS_2") + .setP0(35) + .setQ0(20) + .add(); + + VoltageLevel vl3 = substation.newVoltageLevel() + .setId("VL_3") + .setNominalV(132.0) + .setLowVoltageLimit(118.8) + .setHighVoltageLimit(145.2) + .setTopologyKind(TopologyKind.BUS_BREAKER) + .add(); + vl3.getBusBreakerView().newBus() + .setId("BUS_3") + .add(); + Generator g3 = vl3.newGenerator() + .setId("GEN_3") + .setBus("BUS_3") + .setMinP(0.0) + .setMaxP(140) + .setTargetP(15) + .setTargetV(130) + .setVoltageRegulatorOn(true) + .add(); + g3.newMinMaxReactiveLimits().setMinQ(-3000.0).setMaxQ(3000.0).add(); + + network.newLine() + .setId("LINE_12") + .setBus1("BUS_1") + .setBus2("BUS_2") + .setR(1.05) + .setX(10.0) + .setG1(0.0000005) + .add(); + network.newLine() + .setId("LINE_23") + .setBus1("BUS_2") + .setBus2("BUS_3") + .setR(1.05) + .setX(10.0) + .setG1(0.0000005) + .add(); + + OpenLoadFlowParameters.get(parameters) +// .setVoltageInitModeOverride(OpenLoadFlowParameters.VoltageInitModeOverride.FULL_VOLTAGE) + .setSlackBusSelectionMode(SlackBusSelectionMode.NAME) + .setVoltageRemoteControl(false) + .setSlackBusId("VL_1_0"); + + // Line Characteristics in per-unit + double rPU = 0.0; + double xPU = 1e-5; + // Voltage Mismatch + double alpha = 0.75; + PerturbationFactory.VoltagePerturbation perturbation = PerturbationFactory.getVoltagePerturbation(network); + System.out.println(perturbation); + PerturbationFactory.applyVoltagePerturbation(network, perturbation, rPU, xPU, alpha); + + fixReacLim(network, listMinQ, listMaxQ); +// OpenLoadFlowParameters.get(parameters).setAcSolverType("NEWTON_RAPHSON"); + LoadFlowResult result = loadFlowRunner.run(network, parameters); + assertTrue(result.isFullyConverged(), "Not Fully Converged"); + ReacLimitsTestsUtils.checkSwitches(network, listMinQ, listMaxQ); + ReacLimitsTestsUtils.verifNewtonRaphson(network, parameters, loadFlowRunner, 20); + } + + @Test + public void ieee14ActivePowerPerturbed() { + String perturbation = "ActivePowerLocal"; + logFile = "D:\\Documents\\Logs_Tests\\Logs_ieee14_" + perturbation + ".txt"; + Network network = IeeeCdfNetworkFactory.create14(); + testprocess(logFile, network, perturbation, 1.2); + ReacLimitsTestsUtils.verifNewtonRaphson(network, parameters, loadFlowRunner, 20); + } + + @Test + public void ieee30ActivePowerPerturbed() { + String perturbation = "ActivePowerLocal"; + logFile = "D:\\Documents\\Logs_Tests\\Logs_ieee30_" + perturbation + ".txt"; + Network network = IeeeCdfNetworkFactory.create30(); + testprocess(logFile, network, perturbation, 1.2); + ReacLimitsTestsUtils.verifNewtonRaphson(network, parameters, loadFlowRunner, 20); + } + + @Test + public void ieee30ReactivePowerPerturbed() { + String perturbation = "ReactivePower"; + logFile = "D:\\Documents\\Logs_Tests\\Logs_ieee30_" + perturbation + ".txt"; + Network network = IeeeCdfNetworkFactory.create30(); + testprocess(logFile, network, perturbation, 1E11); + ReacLimitsTestsUtils.verifNewtonRaphson(network, parameters, loadFlowRunner, 20); + } + + @Test + public void ieee30VoltagePerturbed() { + // set up network + HashMap listMinQ = new HashMap<>(); + HashMap listMaxQ = new HashMap<>(); + parameters.setUseReactiveLimits(true); + Network network = IeeeCdfNetworkFactory.create30(); + fixReacLim(network, listMinQ, listMaxQ); + + // Line Characteristics in per-unit + double rPU = 0.0; + double xPU = 1e-5; + // Voltage Mismatch + double alpha = 0.95; + PerturbationFactory.VoltagePerturbation perturbation = PerturbationFactory.getVoltagePerturbation(network); + PerturbationFactory.applyVoltagePerturbation(network, perturbation, rPU, xPU, alpha); + + // solve and check + LoadFlowResult result = loadFlowRunner.run(network, parameters); + assertTrue(result.isFullyConverged(), "Not Fully Converged"); + ReacLimitsTestsUtils.checkSwitches(network, listMinQ, listMaxQ); + ReacLimitsTestsUtils.verifNewtonRaphson(network, parameters, loadFlowRunner, 20); + } + + @Test + public void ieee118ActivePowerPerturbed() { + String perturbation = "ActivePowerLocal"; + logFile = "D:\\Documents\\Logs_Tests\\Logs_ieee118_" + perturbation + ".txt"; + Network network = IeeeCdfNetworkFactory.create118(); + testprocess(logFile, network, perturbation, 1.2); + ReacLimitsTestsUtils.verifNewtonRaphson(network, parameters, loadFlowRunner, 20); + } + + @Test + public void ieee118ReactivePowerPerturbed() { + String perturbation = "ReactivePower"; + logFile = "D:\\Documents\\Logs_Tests\\Logs_ieee118_" + perturbation + ".txt"; + Network network = IeeeCdfNetworkFactory.create118(); + testprocess(logFile, network, perturbation, 1E10); + ReacLimitsTestsUtils.verifNewtonRaphson(network, parameters, loadFlowRunner, 20); + } + + @Test + public void ieee118VoltagePerturbed() { + // set up network + HashMap listMinQ = new HashMap<>(); + HashMap listMaxQ = new HashMap<>(); + parameters.setUseReactiveLimits(true); + Network network = IeeeCdfNetworkFactory.create118(); + fixReacLim(network, listMinQ, listMaxQ); + + // Line Characteristics in per-unit + double rPU = 0.0; + double xPU = 1e-5; + // Voltage Mismatch + double alpha = 0.95; + PerturbationFactory.VoltagePerturbation perturbation = PerturbationFactory.getVoltagePerturbation(network); + PerturbationFactory.applyVoltagePerturbation(network, perturbation, rPU, xPU, alpha); + + // solve and check + LoadFlowResult result = loadFlowRunner.run(network, parameters); + assertTrue(result.isFullyConverged(), "Not Fully Converged"); + ReacLimitsTestsUtils.checkSwitches(network, listMinQ, listMaxQ); + ReacLimitsTestsUtils.verifNewtonRaphson(network, parameters, loadFlowRunner, 20); + } + + @Test + public void ieee300ActivePowerPerturbed() { + String perturbation = "ActivePowerGlobal"; + logFile = "D:\\Documents\\Logs_Tests\\Logs_ieee300_" + perturbation + ".txt"; + Network network = IeeeCdfNetworkFactory.create300(); + testprocess(logFile, network, perturbation, 1.2); + ReacLimitsTestsUtils.verifNewtonRaphson(network, parameters, loadFlowRunner, 20); + } + + @Test + public void ieee300ReactivePowerPerturbed() { + String perturbation = "ReactivePower"; + logFile = "D:\\Documents\\Logs_Tests\\Logs_ieee300_" + perturbation + ".txt"; + Network network = IeeeCdfNetworkFactory.create300(); + testprocess(logFile, network, perturbation, 1E11); + ReacLimitsTestsUtils.verifNewtonRaphson(network, parameters, loadFlowRunner, 20); + } + + @Test + void testxiidm1888ActivePowerOneLoadPerturbed() throws IOException { + String perturbation = "ActivePowerLocal"; + logFile = "D:\\Documents\\Logs_Tests\\Logs_Rte1888_" + perturbation + ".txt"; + Network network = Network.read("D:\\Documents\\Réseaux\\rte1888.xiidm"); + testprocess(logFile, network, perturbation, 1.2); + ReacLimitsTestsUtils.verifNewtonRaphson(network, parameters, loadFlowRunner, 20); + } + + @Test + void testxiidm6515() throws IOException { + String perturbation = "None"; + logFile = "D:\\Documents\\Logs_Tests\\Logs_Rte6515_" + perturbation + ".txt"; + Network network = Network.read("D:\\Documents\\Réseaux\\rte6515.xiidm"); + testprocess(logFile, network, perturbation, 1.2); + } +} diff --git a/src/test/java/com/powsybl/openloadflow/knitro/solver/ReactiveNoJacobienneTest.java b/src/test/java/com/powsybl/openloadflow/knitro/solver/ReactiveNoJacobienneTest.java new file mode 100644 index 0000000..064acdb --- /dev/null +++ b/src/test/java/com/powsybl/openloadflow/knitro/solver/ReactiveNoJacobienneTest.java @@ -0,0 +1,438 @@ +/** + * 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.*; + +import static com.powsybl.openloadflow.util.LoadFlowAssert.assertReactivePowerEquals; +import static org.junit.jupiter.api.Assertions.*; + +/** + * @author Pierre Arvy {@literal } + * @author Yoann Anezin {@literal } + */ +public class ReactiveNoJacobienneTest { + private static final double DEFAULT_TOLERANCE = 1e-2; + private LoadFlow.Runner loadFlowRunner; + private LoadFlowParameters parameters; + + @BeforeEach + void setUp() { + loadFlowRunner = new LoadFlow.Runner(new OpenLoadFlowProvider(new DenseMatrixFactory())); + parameters = new LoadFlowParameters().setUseReactiveLimits(true) + .setDistributedSlack(false); + KnitroLoadFlowParameters knitroLoadFlowParameters = new KnitroLoadFlowParameters(); // set gradient computation mode + knitroLoadFlowParameters.setGradientComputationMode(2); + knitroLoadFlowParameters.setMaxIterations(300); + knitroLoadFlowParameters.setKnitroSolverType(KnitroSolverParameters.KnitroSolverType.REACTIVLIMITS); + parameters.addExtension(KnitroLoadFlowParameters.class, knitroLoadFlowParameters); + //parameters.setVoltageInitMode(LoadFlowParameters.VoltageInitMode.DC_VALUES); + OpenLoadFlowParameters.create(parameters).setAcSolverType(KnitroSolverFactory.NAME); + OpenLoadFlowParameters.get(parameters).setVoltageInitModeOverride(OpenLoadFlowParameters.VoltageInitModeOverride.FULL_VOLTAGE); + } + + @Test + void testReacLimEurostagQlow() { + HashMap listMinQ = new HashMap<>(); + HashMap listMaxQ = new HashMap<>(); + Network network = EurostagFactory.fix(EurostagTutorialExample1Factory.create()); + + // access to already created equipments + Load load = network.getLoad("LOAD"); + VoltageLevel vlgen = network.getVoltageLevel("VLGEN"); + TwoWindingsTransformer nhv2Nload = network.getTwoWindingsTransformer("NHV2_NLOAD"); + Generator gen = network.getGenerator("GEN"); + Substation p1 = network.getSubstation("P1"); + + // reduce GEN reactive range + gen.newMinMaxReactiveLimits() + .setMinQ(0) + .setMaxQ(280) + .add(); + listMinQ.put(gen.getId(), -280.0); + listMaxQ.put(gen.getId(), 280.0); + + // 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(); + Generator gen2 = vlgen2.newGenerator() + .setId("GEN2") + .setBus("NGEN2") + .setConnectableBus("NGEN2") + .setMinP(-9999.99) + .setMaxP(9999.99) + .setVoltageRegulatorOn(true) + .setTargetV(24.5) + .setTargetP(100) + .add(); + gen2.newMinMaxReactiveLimits() + .setMinQ(250) + .setMaxQ(300) + .add(); + listMinQ.put(gen2.getId(), 250.0); + listMaxQ.put(gen2.getId(), 300.0); + int zb380 = 380 * 380 / 100; + TwoWindingsTransformer ngen2Nhv1 = p1.newTwoWindingsTransformer() + .setId("NGEN2_NHV1") + .setBus1("NGEN2") + .setConnectableBus1("NGEN2") + .setRatedU1(24.0) + .setBus2("NHV1") + .setConnectableBus2("NHV1") + .setRatedU2(400.0) + .setR(0.24 / 1800 * zb380) + .setX(Math.sqrt(10 * 10 - 0.24 * 0.24) / 1800 * zb380) + .add(); + + // fix active power balance + load.setP0(699.838); + + LoadFlowResult result = loadFlowRunner.run(network, parameters); + assertTrue(result.isFullyConverged()); +// assertReactivePowerEquals(-8.094, gen.getTerminal()); + assertReactivePowerEquals(-250, gen2.getTerminal()); // GEN is correctly limited to 250 MVar + assertReactivePowerEquals(250, ngen2Nhv1.getTerminal1()); + assertReactivePowerEquals(-200, nhv2Nload.getTerminal2()); + ReacLimitsTestsUtils.checkSwitches(network, listMinQ, listMaxQ); + ReacLimitsTestsUtils.verifNewtonRaphson(network, parameters, loadFlowRunner, 0); + } + + @Test + void testReacLimEurostagQup() { + HashMap listMinQ = new HashMap<>(); + HashMap listMaxQ = new HashMap<>(); + Network network = EurostagFactory.fix(EurostagTutorialExample1Factory.create()); + + // access to already created equipments + Load load = network.getLoad("LOAD"); + + VoltageLevel vlgen = network.getVoltageLevel("VLGEN"); + TwoWindingsTransformer nhv2Nload = network.getTwoWindingsTransformer("NHV2_NLOAD"); + Generator gen = network.getGenerator("GEN"); + Substation p1 = network.getSubstation("P1"); + + // reduce GEN reactive range + gen.newMinMaxReactiveLimits() + .setMinQ(0) + .setMaxQ(280) + .add(); + listMinQ.put(gen.getId(), -280.0); + listMaxQ.put(gen.getId(), 280.0); + + // 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(); + Generator gen2 = vlgen2.newGenerator() + .setId("GEN2") + .setBus("NGEN2") + .setConnectableBus("NGEN2") + .setMinP(-9999.99) + .setMaxP(9999.99) + .setVoltageRegulatorOn(true) + .setTargetV(24.5) + .setTargetP(100) + .add(); + gen2.newMinMaxReactiveLimits() + .setMinQ(0) + .setMaxQ(100) + .add(); + listMinQ.put(gen2.getId(), 0.0); + listMaxQ.put(gen2.getId(), 100.0); + int zb380 = 380 * 380 / 100; + TwoWindingsTransformer ngen2Nhv1 = p1.newTwoWindingsTransformer() + .setId("NGEN2_NHV1") + .setBus1("NGEN2") + .setConnectableBus1("NGEN2") + .setRatedU1(24.0) + .setBus2("NHV1") + .setConnectableBus2("NHV1") + .setRatedU2(400.0) + .setR(0.24 / 1800 * zb380) + .setX(Math.sqrt(10 * 10 - 0.24 * 0.24) / 1800 * zb380) + .add(); + + // fix active power balance + load.setP0(699.838); + + LoadFlowResult result = loadFlowRunner.run(network, parameters); + assertTrue(result.isFullyConverged()); +// assertReactivePowerEquals(-164.315, gen.getTerminal()); +// assertReactivePowerEquals(-100, gen2.getTerminal()); // GEN is correctly limited to 100 MVar +// assertReactivePowerEquals(100, ngen2Nhv1.getTerminal1()); +// assertReactivePowerEquals(-200, nhv2Nload.getTerminal2()); + ReacLimitsTestsUtils.checkSwitches(network, listMinQ, listMaxQ); + ReacLimitsTestsUtils.verifNewtonRaphson(network, parameters, loadFlowRunner, 0); + } + + @Test + void testReacLimEurostagQupWithLoad() { + HashMap listMinQ = new HashMap<>(); + HashMap listMaxQ = new HashMap<>(); + Network network = EurostagFactory.fix(EurostagTutorialExample1Factory.create()); + + // access to already created equipments + Load load = network.getLoad("LOAD"); + + VoltageLevel vlgen = network.getVoltageLevel("VLGEN"); + TwoWindingsTransformer nhv2Nload = network.getTwoWindingsTransformer("NHV2_NLOAD"); + Generator gen = network.getGenerator("GEN"); + Substation p1 = network.getSubstation("P1"); + + // reduce GEN reactive range + gen.newMinMaxReactiveLimits() + .setMinQ(0) + .setMaxQ(280) + .add(); + listMinQ.put(gen.getId(), -280.0); + listMaxQ.put(gen.getId(), 280.0); + + // 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(); + Generator gen2 = vlgen2.newGenerator() + .setId("GEN2") + .setBus("NGEN2") + .setConnectableBus("NGEN2") + .setMinP(-9999.99) + .setMaxP(9999.99) + .setVoltageRegulatorOn(true) + .setTargetV(24.5) + .setTargetP(100) + .add(); + gen2.newMinMaxReactiveLimits() + .setMinQ(0) + .setMaxQ(100) + .add(); + listMinQ.put(gen2.getId(), 0.0); + listMaxQ.put(gen2.getId(), 100.0); + Load load2 = vlgen2.newLoad() + .setId("LOAD2") + .setBus("NGEN2") + .setConnectableBus("NGEN2") + .setP0(0.0) + .setQ0(30.0) + .add(); + int zb380 = 380 * 380 / 100; + TwoWindingsTransformer ngen2Nhv1 = p1.newTwoWindingsTransformer() + .setId("NGEN2_NHV1") + .setBus1("NGEN2") + .setConnectableBus1("NGEN2") + .setRatedU1(24.0) + .setBus2("NHV1") + .setConnectableBus2("NHV1") + .setRatedU2(400.0) + .setR(0.24 / 1800 * zb380) + .setX(Math.sqrt(10 * 10 - 0.24 * 0.24) / 1800 * zb380) + .add(); + + // fix active power balance + load.setP0(699.838); + + LoadFlowResult result = loadFlowRunner.run(network, parameters); + assertTrue(result.isFullyConverged()); + //assertReactivePowerEquals(-196.263, gen.getTerminal()); + assertReactivePowerEquals(-100, gen2.getTerminal()); // GEN is correctly limited to 100 MVar + assertReactivePowerEquals(70, ngen2Nhv1.getTerminal1()); + assertReactivePowerEquals(-200, nhv2Nload.getTerminal2()); + ReacLimitsTestsUtils.checkSwitches(network, listMinQ, listMaxQ); + ReacLimitsTestsUtils.verifNewtonRaphson(network, parameters, loadFlowRunner, 0); + } + + @Test + void testReacLimEurostagQupWithGen() { + HashMap listMinQ = new HashMap<>(); + HashMap listMaxQ = new HashMap<>(); + Network network = EurostagFactory.fix(EurostagTutorialExample1Factory.create()); + + // access to already created equipments + Load load = network.getLoad("LOAD"); + + VoltageLevel vlgen = network.getVoltageLevel("VLGEN"); + TwoWindingsTransformer nhv2Nload = network.getTwoWindingsTransformer("NHV2_NLOAD"); + Generator gen = network.getGenerator("GEN"); + Substation p1 = network.getSubstation("P1"); + + // reduce GEN reactive range + gen.newMinMaxReactiveLimits() + .setMinQ(0) + .setMaxQ(280) + .add(); + listMinQ.put(gen.getId(), -280.0); + listMaxQ.put(gen.getId(), 280.0); + + // 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(); + Generator gen2 = vlgen2.newGenerator() + .setId("GEN2") + .setBus("NGEN2") + .setConnectableBus("NGEN2") + .setMinP(-9999.99) + .setMaxP(9999.99) + .setVoltageRegulatorOn(true) + .setTargetV(24.5) + .setTargetP(100) + .add(); + gen2.newMinMaxReactiveLimits() + .setMinQ(0) + .setMaxQ(100) + .add(); + listMinQ.put(gen2.getId(), 0.0); + listMaxQ.put(gen2.getId(), 100.0); + Generator gen2Bis = vlgen2.newGenerator() + .setId("GEN2BIS") + .setBus("NGEN2") + .setConnectableBus("NGEN2") + .setMinP(-9999.99) + .setMaxP(9999.99) + .setVoltageRegulatorOn(true) + .setTargetV(24.5) + .setTargetP(50) + .add(); + gen2Bis.newMinMaxReactiveLimits() + .setMinQ(0) + .setMaxQ(40) + .add(); + listMinQ.put(gen2Bis.getId(), 0.0); + listMaxQ.put(gen2Bis.getId(), 40.0); + int zb380 = 380 * 380 / 100; + TwoWindingsTransformer ngen2Nhv1 = p1.newTwoWindingsTransformer() + .setId("NGEN2_NHV1") + .setBus1("NGEN2") + .setConnectableBus1("NGEN2") + .setRatedU1(24.0) + .setBus2("NHV1") + .setConnectableBus2("NHV1") + .setRatedU2(400.0) + .setR(0.24 / 1800 * zb380) + .setX(Math.sqrt(10 * 10 - 0.24 * 0.24) / 1800 * zb380) + .add(); + + // fix active power balance + load.setP0(699.838); + LoadFlowResult 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()); + ReacLimitsTestsUtils.checkSwitches(network, listMinQ, listMaxQ); + ReacLimitsTestsUtils.verifNewtonRaphson(network, parameters, loadFlowRunner, 0); + } + + @Test + void testReacLimIeee14() { + HashMap listMinQ = new HashMap<>(); + HashMap listMaxQ = new HashMap<>(); + parameters.setUseReactiveLimits(true); + Network network = IeeeCdfNetworkFactory.create14(); + for (var g : network.getGenerators()) { + if (g.getReactiveLimits().getMinQ(g.getTerminal().getBusView().getBus().getP()) > -1.7976931348623157E308) { + listMinQ.put(g.getId(), g.getReactiveLimits().getMinQ(g.getTerminal().getBusView().getBus().getP())); + listMaxQ.put(g.getId(), g.getReactiveLimits().getMaxQ(g.getTerminal().getBusView().getBus().getP())); + } + } + LoadFlowResult result = loadFlowRunner.run(network, parameters); + assertTrue(result.isFullyConverged(), "Not Fully Converged"); + ReacLimitsTestsUtils.checkSwitches(network, listMinQ, listMaxQ); + ReacLimitsTestsUtils.verifNewtonRaphson(network, parameters, loadFlowRunner, 0); + } + + @Test + void testReacLimIeee30() { + HashMap listMinQ = new HashMap<>(); + HashMap listMaxQ = new HashMap<>(); + parameters.setUseReactiveLimits(true); + Network network = IeeeCdfNetworkFactory.create30(); + for (var g : network.getGenerators()) { + if (g.getReactiveLimits().getMinQ(g.getTerminal().getBusView().getBus().getP()) > -1.7976931348623157E308) { + listMinQ.put(g.getId(), g.getReactiveLimits().getMinQ(g.getTerminal().getBusView().getBus().getP())); + listMaxQ.put(g.getId(), g.getReactiveLimits().getMaxQ(g.getTerminal().getBusView().getBus().getP())); + } + } + LoadFlowResult result = loadFlowRunner.run(network, parameters); + assertTrue(result.isFullyConverged(), "Not Fully Converged"); + ReacLimitsTestsUtils.checkSwitches(network, listMinQ, listMaxQ); + ReacLimitsTestsUtils.verifNewtonRaphson(network, parameters, loadFlowRunner, 0); + } + + @Test + void testReacLimIeee118() { + HashMap listMinQ = new HashMap<>(); + HashMap listMaxQ = new HashMap<>(); + parameters.setUseReactiveLimits(true); + Network network = IeeeCdfNetworkFactory.create118(); + for (var g : network.getGenerators()) { + if (g.getReactiveLimits().getMinQ(g.getTerminal().getBusView().getBus().getP()) > -1.7976931348623157E308) { + listMinQ.put(g.getId(), g.getReactiveLimits().getMinQ(g.getTerminal().getBusView().getBus().getP())); + listMaxQ.put(g.getId(), g.getReactiveLimits().getMaxQ(g.getTerminal().getBusView().getBus().getP())); + } + } + LoadFlowResult result = loadFlowRunner.run(network, parameters); + assertTrue(result.isFullyConverged(), "Not Fully Converged"); + ReacLimitsTestsUtils.checkSwitches(network, listMinQ, listMaxQ); + ReacLimitsTestsUtils.verifNewtonRaphson(network, parameters, loadFlowRunner, 0); + } + + @Test + void testReacLimIeee300() { + HashMap listMinQ = new HashMap<>(); + HashMap listMaxQ = new HashMap<>(); + parameters.setUseReactiveLimits(true); + Network network = IeeeCdfNetworkFactory.create300(); + for (var g : network.getGenerators()) { + if (g.getReactiveLimits().getMinQ(g.getTerminal().getBusView().getBus().getP()) > -1.7976931348623157E308) { + listMinQ.put(g.getId(), g.getReactiveLimits().getMinQ(g.getTerminal().getBusView().getBus().getP())); + listMaxQ.put(g.getId(), g.getReactiveLimits().getMaxQ(g.getTerminal().getBusView().getBus().getP())); + } + } + LoadFlowResult result = loadFlowRunner.run(network, parameters); + assertTrue(result.isFullyConverged(), "Not Fully Converged"); + ReacLimitsTestsUtils.checkSwitches(network, listMinQ, listMaxQ); + ReacLimitsTestsUtils.verifNewtonRaphson(network, parameters, loadFlowRunner, 0); + } +} diff --git a/src/test/java/com/powsybl/openloadflow/knitro/solver/ReactiveWithJacobienneTest.java b/src/test/java/com/powsybl/openloadflow/knitro/solver/ReactiveWithJacobienneTest.java new file mode 100644 index 0000000..1e4cea0 --- /dev/null +++ b/src/test/java/com/powsybl/openloadflow/knitro/solver/ReactiveWithJacobienneTest.java @@ -0,0 +1,456 @@ +/** + * 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.SparseMatrixFactory; +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.*; + +import static com.powsybl.openloadflow.util.LoadFlowAssert.assertReactivePowerEquals; +import static org.junit.jupiter.api.Assertions.*; + +/** + * @author Pierre Arvy {@literal } + * @author Yoann Anezin {@literal } + */ +public class ReactiveWithJacobienneTest { + private static final double DEFAULT_TOLERANCE = 1e-3; + private LoadFlow.Runner loadFlowRunner; + private LoadFlowParameters parameters; + + @BeforeEach + void setUp() { + loadFlowRunner = new LoadFlow.Runner(new OpenLoadFlowProvider(new SparseMatrixFactory())); + parameters = new LoadFlowParameters().setUseReactiveLimits(true) + .setDistributedSlack(false); + KnitroLoadFlowParameters knitroLoadFlowParameters = new KnitroLoadFlowParameters(); // set gradient computation mode + knitroLoadFlowParameters.setGradientComputationMode(1); + knitroLoadFlowParameters.setMaxIterations(300); + knitroLoadFlowParameters.setKnitroSolverType(KnitroSolverParameters.KnitroSolverType.REACTIVLIMITS); + parameters.addExtension(KnitroLoadFlowParameters.class, knitroLoadFlowParameters); + //parameters.setVoltageInitMode(LoadFlowParameters.VoltageInitMode.DC_VALUES); + //OpenLoadFlowParameters.create(parameters).setAcSolverType("NEWTON_RAPHSON"); + OpenLoadFlowParameters.create(parameters).setAcSolverType(KnitroSolverFactory.NAME); + OpenLoadFlowParameters.get(parameters).setVoltageInitModeOverride(OpenLoadFlowParameters.VoltageInitModeOverride.FULL_VOLTAGE); + } + + @Test + void testReacLimEurostagQlow() { + HashMap listMinQ = new HashMap<>(); + HashMap listMaxQ = new HashMap<>(); + Network network = EurostagFactory.fix(EurostagTutorialExample1Factory.create()); + + // access to already created equipments + Load load = network.getLoad("LOAD"); + VoltageLevel vlgen = network.getVoltageLevel("VLGEN"); + TwoWindingsTransformer nhv2Nload = network.getTwoWindingsTransformer("NHV2_NLOAD"); + Generator gen = network.getGenerator("GEN"); + Substation p1 = network.getSubstation("P1"); + + // reduce GEN reactive range + gen.newMinMaxReactiveLimits() + .setMinQ(0) + .setMaxQ(280) + .add(); + listMinQ.put(gen.getId(), -280.0); + listMaxQ.put(gen.getId(), 280.0); + + // 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(); + Generator gen2 = vlgen2.newGenerator() + .setId("GEN2") + .setBus("NGEN2") + .setConnectableBus("NGEN2") + .setMinP(-9999.99) + .setMaxP(9999.99) + .setVoltageRegulatorOn(true) + .setTargetV(24.5) + .setTargetP(100) + .add(); + gen2.newMinMaxReactiveLimits() + .setMinQ(250) + .setMaxQ(300) + .add(); + listMinQ.put(gen2.getId(), 250.0); + listMaxQ.put(gen2.getId(), 300.0); + int zb380 = 380 * 380 / 100; + TwoWindingsTransformer ngen2Nhv1 = p1.newTwoWindingsTransformer() + .setId("NGEN2_NHV1") + .setBus1("NGEN2") + .setConnectableBus1("NGEN2") + .setRatedU1(24.0) + .setBus2("NHV1") + .setConnectableBus2("NHV1") + .setRatedU2(400.0) + .setR(0.24 / 1800 * zb380) + .setX(Math.sqrt(10 * 10 - 0.24 * 0.24) / 1800 * zb380) + .add(); + + // fix active power balance + load.setP0(699.838); + + LoadFlowResult result = loadFlowRunner.run(network, parameters); + assertTrue(result.isFullyConverged()); + //assertReactivePowerEquals(-8.094, gen.getTerminal()); + assertReactivePowerEquals(-250, gen2.getTerminal()); // GEN is correctly limited to 250 MVar + assertReactivePowerEquals(250, ngen2Nhv1.getTerminal1()); + assertReactivePowerEquals(-200, nhv2Nload.getTerminal2()); + ReacLimitsTestsUtils.checkSwitches(network, listMinQ, listMaxQ); + ReacLimitsTestsUtils.verifNewtonRaphson(network, parameters, loadFlowRunner, 0); + } + + @Test + void testReacLimEurostagQup() { + HashMap listMinQ = new HashMap<>(); + HashMap listMaxQ = new HashMap<>(); + Network network = EurostagFactory.fix(EurostagTutorialExample1Factory.create()); + + // access to already created equipments + Load load = network.getLoad("LOAD"); + VoltageLevel vlgen = network.getVoltageLevel("VLGEN"); + TwoWindingsTransformer nhv2Nload = network.getTwoWindingsTransformer("NHV2_NLOAD"); + Generator gen = network.getGenerator("GEN"); + Substation p1 = network.getSubstation("P1"); + + // reduce GEN reactive range + gen.newMinMaxReactiveLimits() + .setMinQ(0) + .setMaxQ(280) + .add(); + listMinQ.put(gen.getId(), -280.0); + listMaxQ.put(gen.getId(), 280.0); + + // 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(); + Generator gen2 = vlgen2.newGenerator() + .setId("GEN2") + .setBus("NGEN2") + .setConnectableBus("NGEN2") + .setMinP(-9999.99) + .setMaxP(9999.99) + .setVoltageRegulatorOn(true) + .setTargetV(24.5) + .setTargetP(100) + .add(); + gen2.newMinMaxReactiveLimits() + .setMinQ(0) + .setMaxQ(100) + .add(); + listMinQ.put(gen2.getId(), 0.0); + listMaxQ.put(gen2.getId(), 100.0); + int zb380 = 380 * 380 / 100; + TwoWindingsTransformer ngen2Nhv1 = p1.newTwoWindingsTransformer() + .setId("NGEN2_NHV1") + .setBus1("NGEN2") + .setConnectableBus1("NGEN2") + .setRatedU1(24.0) + .setBus2("NHV1") + .setConnectableBus2("NHV1") + .setRatedU2(400.0) + .setR(0.24 / 1800 * zb380) + .setX(Math.sqrt(10 * 10 - 0.24 * 0.24) / 1800 * zb380) + .add(); + + // fix active power balance + load.setP0(699.838); + + parameters.setUseReactiveLimits(true); + LoadFlowResult result = loadFlowRunner.run(network, parameters); + assertTrue(result.isFullyConverged()); + assertReactivePowerEquals(-164.3169, gen.getTerminal()); + assertReactivePowerEquals(-100, gen2.getTerminal()); // GEN is correctly limited to 100 MVar + assertReactivePowerEquals(100, ngen2Nhv1.getTerminal1()); + assertReactivePowerEquals(-200, nhv2Nload.getTerminal2()); + ReacLimitsTestsUtils.checkSwitches(network, listMinQ, listMaxQ); + ReacLimitsTestsUtils.verifNewtonRaphson(network, parameters, loadFlowRunner, 0); + } + + @Test + void testReacLimEurostagQupWithLoad() { + HashMap listMinQ = new HashMap<>(); + HashMap listMaxQ = new HashMap<>(); + Network network = EurostagFactory.fix(EurostagTutorialExample1Factory.create()); + + // access to already created equipments + Load load = network.getLoad("LOAD"); + + VoltageLevel vlgen = network.getVoltageLevel("VLGEN"); + TwoWindingsTransformer nhv2Nload = network.getTwoWindingsTransformer("NHV2_NLOAD"); + Generator gen = network.getGenerator("GEN"); + Substation p1 = network.getSubstation("P1"); + + // reduce GEN reactive range + gen.newMinMaxReactiveLimits() + .setMinQ(0) + .setMaxQ(280) + .add(); + listMinQ.put(gen.getId(), -280.0); + listMaxQ.put(gen.getId(), 280.0); + + // 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(); + Generator gen2 = vlgen2.newGenerator() + .setId("GEN2") + .setBus("NGEN2") + .setConnectableBus("NGEN2") + .setMinP(-9999.99) + .setMaxP(9999.99) + .setVoltageRegulatorOn(true) + .setTargetV(24.5) + .setTargetP(100) + .add(); + gen2.newMinMaxReactiveLimits() + .setMinQ(0) + .setMaxQ(100) + .add(); + listMinQ.put(gen2.getId(), 0.0); + listMaxQ.put(gen2.getId(), 100.0); + Load load2 = vlgen2.newLoad() + .setId("LOAD2") + .setBus("NGEN2") + .setConnectableBus("NGEN2") + .setP0(0.0) + .setQ0(30.0) + .add(); + int zb380 = 380 * 380 / 100; + TwoWindingsTransformer ngen2Nhv1 = p1.newTwoWindingsTransformer() + .setId("NGEN2_NHV1") + .setBus1("NGEN2") + .setConnectableBus1("NGEN2") + .setRatedU1(24.0) + .setBus2("NHV1") + .setConnectableBus2("NHV1") + .setRatedU2(400.0) + .setR(0.24 / 1800 * zb380) + .setX(Math.sqrt(10 * 10 - 0.24 * 0.24) / 1800 * zb380) + .add(); + + // fix active power balance + load.setP0(699.838); + + LoadFlowResult result = loadFlowRunner.run(network, parameters); + assertTrue(result.isFullyConverged()); + assertReactivePowerEquals(-196.264, gen.getTerminal()); + assertReactivePowerEquals(-100, gen2.getTerminal()); // GEN is correctly limited to 100 MVar + assertReactivePowerEquals(70, ngen2Nhv1.getTerminal1()); + assertReactivePowerEquals(-200, nhv2Nload.getTerminal2()); + ReacLimitsTestsUtils.checkSwitches(network, listMinQ, listMaxQ); + ReacLimitsTestsUtils.verifNewtonRaphson(network, parameters, loadFlowRunner, 0); + } + + @Test + void testReacLimEurostagQupWithGen() { + HashMap listMinQ = new HashMap<>(); + HashMap listMaxQ = new HashMap<>(); + Network network = EurostagFactory.fix(EurostagTutorialExample1Factory.create()); + + // access to already created equipments + Load load = network.getLoad("LOAD"); + + VoltageLevel vlgen = network.getVoltageLevel("VLGEN"); + TwoWindingsTransformer nhv2Nload = network.getTwoWindingsTransformer("NHV2_NLOAD"); + Generator gen = network.getGenerator("GEN"); + Substation p1 = network.getSubstation("P1"); + + // reduce GEN reactive range + gen.newMinMaxReactiveLimits() + .setMinQ(0) + .setMaxQ(280) + .add(); + listMinQ.put(gen.getId(), -280.0); + listMaxQ.put(gen.getId(), 280.0); + + // 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(); + Generator gen2 = vlgen2.newGenerator() + .setId("GEN2") + .setBus("NGEN2") + .setConnectableBus("NGEN2") + .setMinP(-9999.99) + .setMaxP(9999.99) + .setVoltageRegulatorOn(true) + .setTargetV(24.5) + .setTargetP(100) + .add(); + gen2.newMinMaxReactiveLimits() + .setMinQ(0) + .setMaxQ(100) + .add(); + listMinQ.put(gen2.getId(), 0.0); + listMaxQ.put(gen2.getId(), 100.0); + Generator gen2Bis = vlgen2.newGenerator() + .setId("GEN2BIS") + .setBus("NGEN2") + .setConnectableBus("NGEN2") + .setMinP(-9999.99) + .setMaxP(9999.99) + .setVoltageRegulatorOn(true) + .setTargetV(24.5) + .setTargetP(50) + .add(); + gen2Bis.newMinMaxReactiveLimits() + .setMinQ(0) + .setMaxQ(40) + .add(); + listMinQ.put(gen2Bis.getId(), 0.0); + listMaxQ.put(gen2Bis.getId(), 40.0); + int zb380 = 380 * 380 / 100; + TwoWindingsTransformer ngen2Nhv1 = p1.newTwoWindingsTransformer() + .setId("NGEN2_NHV1") + .setBus1("NGEN2") + .setConnectableBus1("NGEN2") + .setRatedU1(24.0) + .setBus2("NHV1") + .setConnectableBus2("NHV1") + .setRatedU2(400.0) + .setR(0.24 / 1800 * zb380) + .setX(Math.sqrt(10 * 10 - 0.24 * 0.24) / 1800 * zb380) + .add(); + + // fix active power balance + load.setP0(699.838); + LoadFlowResult result = loadFlowRunner.run(network, parameters); + assertTrue(result.isFullyConverged()); + //assertReactivePowerEquals(-122.715, gen.getTerminal()); + assertReactivePowerEquals(-100, gen2.getTerminal()); // GEN is correctly limited to 100 MVar + assertReactivePowerEquals(140.0, ngen2Nhv1.getTerminal1()); + assertReactivePowerEquals(-200, nhv2Nload.getTerminal2()); + ReacLimitsTestsUtils.checkSwitches(network, listMinQ, listMaxQ); + ReacLimitsTestsUtils.verifNewtonRaphson(network, parameters, loadFlowRunner, 0); + } + + @Test + void testReacLimIeee14() { /* Unfeasible Point */ + HashMap listMinQ = new HashMap<>(); + HashMap listMaxQ = new HashMap<>(); + parameters.setUseReactiveLimits(true); + Network network = IeeeCdfNetworkFactory.create14(); + for (var g : network.getGenerators()) { + if (g.getReactiveLimits().getMinQ(g.getTerminal().getBusView().getBus().getP()) > -1.7976931348623157E308) { + listMinQ.put(g.getId(), g.getReactiveLimits().getMinQ(g.getTerminal().getBusView().getBus().getP())); + listMaxQ.put(g.getId(), g.getReactiveLimits().getMaxQ(g.getTerminal().getBusView().getBus().getP())); + } + } + LoadFlowResult result = loadFlowRunner.run(network, parameters); + assertTrue(result.isFullyConverged(), "Not Fully Converged"); + ReacLimitsTestsUtils.checkSwitches(network, listMinQ, listMaxQ); + ReacLimitsTestsUtils.verifNewtonRaphson(network, parameters, loadFlowRunner, 0); + } + + @Test + void testReacLimIeee30() { + HashMap listMinQ = new HashMap<>(); + HashMap listMaxQ = new HashMap<>(); + parameters.setUseReactiveLimits(true); + Network network = IeeeCdfNetworkFactory.create30(); + for (var g : network.getGenerators()) { + if (g.getReactiveLimits().getMinQ(g.getTerminal().getBusView().getBus().getP()) > -1.7976931348623157E308) { + listMinQ.put(g.getId(), g.getReactiveLimits().getMinQ(g.getTerminal().getBusView().getBus().getP())); + listMaxQ.put(g.getId(), g.getReactiveLimits().getMaxQ(g.getTerminal().getBusView().getBus().getP())); + } + } + LoadFlowResult result = loadFlowRunner.run(network, parameters); + assertTrue(result.isFullyConverged(), "Not Fully Converged"); + ReacLimitsTestsUtils.checkSwitches(network, listMinQ, listMaxQ); + ReacLimitsTestsUtils.verifNewtonRaphson(network, parameters, loadFlowRunner, 0); + } + + @Test + void testReacLimIeee118() { + HashMap listMinQ = new HashMap<>(); + HashMap listMaxQ = new HashMap<>(); + parameters.setUseReactiveLimits(true); + Network network = IeeeCdfNetworkFactory.create118(); + for (var g : network.getGenerators()) { + if (g.getReactiveLimits().getMinQ(g.getTerminal().getBusView().getBus().getP()) > -1.7976931348623157E308) { + listMinQ.put(g.getId(), g.getReactiveLimits().getMinQ(g.getTerminal().getBusView().getBus().getP())); + listMaxQ.put(g.getId(), g.getReactiveLimits().getMaxQ(g.getTerminal().getBusView().getBus().getP())); + } + } + LoadFlowResult result = loadFlowRunner.run(network, parameters); + assertTrue(result.isFullyConverged(), "Not Fully Converged"); + ReacLimitsTestsUtils.checkSwitches(network, listMinQ, listMaxQ); + ReacLimitsTestsUtils.verifNewtonRaphson(network, parameters, loadFlowRunner, 0); + } + + @Test + void testReacLimIeee300() { + HashMap listMinQ = new HashMap<>(); + HashMap listMaxQ = new HashMap<>(); + parameters.setUseReactiveLimits(true); + Network network = IeeeCdfNetworkFactory.create300(); + for (var g : network.getGenerators()) { + if (g.getReactiveLimits().getMinQ(g.getTerminal().getBusView().getBus().getP()) > -1.7976931348623157E308) { + listMinQ.put(g.getId(), g.getReactiveLimits().getMinQ(g.getTerminal().getBusView().getBus().getP())); + listMaxQ.put(g.getId(), g.getReactiveLimits().getMaxQ(g.getTerminal().getBusView().getBus().getP())); + } + } + LoadFlowResult result = loadFlowRunner.run(network, parameters); + assertTrue(result.isFullyConverged(), "Not Fully Converged"); + ReacLimitsTestsUtils.checkSwitches(network, listMinQ, listMaxQ); + ReacLimitsTestsUtils.verifNewtonRaphson(network, parameters, loadFlowRunner, 0); + } + + @Test + void testxiidm() { + HashMap listMinQ = new HashMap<>(); + HashMap listMaxQ = new HashMap<>(); + parameters.setUseReactiveLimits(true); + Network network = Network.read("D:\\Documents\\Réseaux\\rte1888.xiidm"); + LoadFlowResult result = loadFlowRunner.run(network, parameters); + assertTrue(result.isFullyConverged(), "Not Fully Converged"); + for (var g : network.getGenerators()) { + if (g.getReactiveLimits().getMinQ(g.getTargetP()) > -1.7976931348623157E308) { + listMinQ.put(g.getId(), g.getReactiveLimits().getMinQ(g.getTargetP())); + listMaxQ.put(g.getId(), g.getReactiveLimits().getMaxQ(g.getTargetP())); + } + } + ReacLimitsTestsUtils.checkSwitches(network, listMinQ, listMaxQ); + ReacLimitsTestsUtils.verifNewtonRaphson(network, parameters, loadFlowRunner, 0); + } +} diff --git a/src/test/java/com/powsybl/openloadflow/knitro/solver/ResilientAcLoadFlowPerturbationTest.java b/src/test/java/com/powsybl/openloadflow/knitro/solver/ResilientAcLoadFlowPerturbationTest.java index 44bd7a4..8ec0ac8 100644 --- a/src/test/java/com/powsybl/openloadflow/knitro/solver/ResilientAcLoadFlowPerturbationTest.java +++ b/src/test/java/com/powsybl/openloadflow/knitro/solver/ResilientAcLoadFlowPerturbationTest.java @@ -54,7 +54,8 @@ void setUp() { private void configureSolver(String solver) { OpenLoadFlowParameters.create(parameters) .setSlackBusSelectionMode(SlackBusSelectionMode.MOST_MESHED) - .setAcSolverType(solver); + .setAcSolverType(solver) + .setVoltageInitModeOverride(OpenLoadFlowParameters.VoltageInitModeOverride.FULL_VOLTAGE); if (RKN.equals(solver)) { KnitroLoadFlowParameters knitroParams = new KnitroLoadFlowParameters(); @@ -124,7 +125,7 @@ void testVoltagePerturbationOnVariousI3ENetworks(NetworkPair pair) { double rPU = 0.0; double xPU = 1e-5; // Voltage Mismatch - double alpha = 0.95; + double alpha = 1.0; voltagePerturbationTest(rknNetwork, nrNetwork, baseFilename, rPU, xPU, alpha); } @@ -204,7 +205,6 @@ void testActivePowerPerturbationOnVariousI3ENetworks(NetworkPair pair) { // Final perturbed load's percentage double alpha = 0.10; - activePowerPerturbationTest(rknNetwork, nrNetwork, baseFilename, alpha); } diff --git a/src/test/java/com/powsybl/openloadflow/knitro/solver/utils/TempoTest.java b/src/test/java/com/powsybl/openloadflow/knitro/solver/utils/TempoTest.java new file mode 100644 index 0000000..73190ae --- /dev/null +++ b/src/test/java/com/powsybl/openloadflow/knitro/solver/utils/TempoTest.java @@ -0,0 +1,101 @@ +package com.powsybl.openloadflow.knitro.solver.utils; + +/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + * This example demonstrates how to use Knitro to solve the following + * simple mathematical program with equilibrium/complementarity + * constraints (MPEC/MPCC). + * + * min (x0 - 5)^2 + (2 x1 + 1)^2 + * s.t. -1.5 x0 + 2 x1 + x2 - 0.5 x3 + x4 = 2 + * x2 complements (3 x0 - x1 - 3) + * x3 complements (-x0 + 0.5 x1 + 4) + * x4 complements (-x0 - x1 + 7) + * x0, x1, x2, x3, x4 >= 0 + * + * The complementarity constraints must be converted so that one + * nonnegative variable complements another nonnegative variable. + * + * min (x0 - 5)^2 + (2 x1 + 1)^2 + * s.t. -1.5 x0 + 2 x1 + x2 - 0.5 x3 + x4 = 2 (c0) + * 3 x0 - x1 - 3 - x5 = 0 (c1) + * -x0 + 0.5 x1 + 4 - x6 = 0 (c2) + * -x0 - x1 + 7 - x7 = 0 (c3) + * x2 complements x5 + * x3 complements x6 + * x4 complements x7 + * x0, x1, x2, x3, x4, x5, x6, x7 >= 0 + * + * The solution is (1, 0, 3.5, 0, 0, 0, 3, 6), with objective value 17. + *++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/ + +import com.artelys.knitro.api.*; + +import java.util.Arrays; + +public final class TempoTest { + + private TempoTest() { + throw new UnsupportedOperationException(); + } + + private static class ProblemMPEC1 extends KNProblem { + public ProblemMPEC1() throws KNException { + super(8, 4); + + // Variables + setVarLoBnds(Arrays.asList(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)); + setXInitial(Arrays.asList(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)); + + // Constraints + setConEqBnds(Arrays.asList(2.0, 3.0, -4.0, -7.0)); + setCompConstraintsTypes(Arrays.asList(KNConstants.KN_CCTYPE_VARVAR, KNConstants.KN_CCTYPE_VARVAR, KNConstants.KN_CCTYPE_VARVAR)); + // x2 x3 and x4 complements respectively x5, x6 and x7 + setCompConstraintsParts(Arrays.asList(2, 3, 4), Arrays.asList(5, 6, 7)); + + // Objective structure + /* Note that the objective (x0 - 5)^2 + (2 x1 + 1)^2 when + * expanded becomes: + * x0^2 + 4 x1^2 - 10 x0 + 4 x1 + 26 */ + /* Add quadratic coefficients for the objective */ + setObjectiveQuadraticPart(Arrays.asList(0, 1), Arrays.asList(0, 1), Arrays.asList(1.0, 4.0)); + /* Add linear coefficients for the objective */ + setObjectiveLinearPart(Arrays.asList(0, 1), Arrays.asList(-10.0, 4.0)); + /* Add constant to the objective */ + setObjConstPart(26.0); + + // Constraints Strucutres + /* c0 */ + addConstraintLinearPart(0, Arrays.asList(0, 1, 2, 3, 4), Arrays.asList(-1.5, 2.0, 1.0, -0.5, 1.0)); + /* c1 */ + addConstraintLinearPart(1, Arrays.asList(0, 1, 5), Arrays.asList(3.0, -1.0, -1.0)); + /* c2 */ + addConstraintLinearPart(2, Arrays.asList(0, 1, 6), Arrays.asList(-1.0, 0.5, -1.0)); + /* c3 */ + addConstraintLinearPart(3, Arrays.asList(0, 1, 7), Arrays.asList(-1.0, -1.0, -1.0)); + } + } + + public static void main(String[] args) throws KNException { + // Create a problem instance. + ProblemMPEC1 instance = new ProblemMPEC1(); + // Create a solver + try (KNSolver solver = new KNSolver(instance)) { + + solver.initProblem(); + solver.solve(); + + KNSolution solution = solver.getSolution(); + + System.out.print("\n\n"); + System.out.format("Knitro converged with final status = %d%n", solution.getStatus()); + System.out.format(" optimal objective value = %f%n", solution.getObjValue()); + System.out.format(" optimal primal values x0=%f%n", solution.getX().get(0)); + System.out.format(" x1=%f%n", solution.getX().get(1)); + System.out.format(" x2=%f complements x5=%f%n", solution.getX().get(2), solution.getX().get(5)); + System.out.format(" x3=%f complements x6=%f%n", solution.getX().get(3), solution.getX().get(6)); + System.out.format(" x4=%f complements x7=%f%n", solution.getX().get(4), solution.getX().get(7)); + System.out.format(" feasibility violation = %f%n", solver.getAbsFeasError()); + System.out.format(" KKT optimality violation = %f%n", solver.getAbsOptError()); + } + } +} diff --git a/src/test/resources/logback-test.xml b/src/test/resources/logback-test.xml index 5416126..9e9e36f 100644 --- a/src/test/resources/logback-test.xml +++ b/src/test/resources/logback-test.xml @@ -13,7 +13,22 @@ %-5p %d{HH:mm:ss.SSS} %-20C{1} | %m%n - + + + D:/Documents/La_Doc/logrte1888.log + + + D:/Documents/La_Doc/logrte1888.%d{yyyy-MM-dd}.log + 30 + + + + %d{yyyy-MM-dd HH:mm:ss} %-5level %logger{36} - %msg%n + + + + +