diff --git a/CHANGELOG.md b/CHANGELOG.md index 41485616d3..af434e6b73 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -25,6 +25,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Updated docker documentation [#488](https://github.com/ie3-institute/simona/issues/488) - Added support classes for transformer tap position calculation [#1543](https://github.com/ie3-institute/simona/issues/1543) - Added basic external em service [#1566](https://github.com/ie3-institute/simona/issues/1566) +- Implement energy limit flex options and adapt optimization [#1572](https://github.com/ie3-institute/simona/issues/1572) ### Changed - Upgraded `scala2` to `scala3` [#53](https://github.com/ie3-institute/simona/issues/53) diff --git a/src/main/scala/edu/ie3/simona/model/em/OptimizedFlexStrat.scala b/src/main/scala/edu/ie3/simona/model/em/OptimizedFlexStrat.scala deleted file mode 100644 index f5f44fe904..0000000000 --- a/src/main/scala/edu/ie3/simona/model/em/OptimizedFlexStrat.scala +++ /dev/null @@ -1,338 +0,0 @@ -/* - * © 2025. TU Dortmund University, - * Institute of Energy Systems, Energy Efficiency and Energy Economics, - * Research group Distribution grid planning and operation - */ - -package edu.ie3.simona.model.em - -import edu.ie3.datamodel.models.input.AssetInput -import edu.ie3.simona.exceptions.CriticalFailureException -import edu.ie3.simona.model.em.OptimizedFlexStrat.* -import edu.ie3.simona.ontology.messages.flex.MathFlexOptions -import edu.ie3.simona.ontology.messages.flex.MathFlexOptions.{ - OperationVars, - SoftConstraint, -} -import optimus.algebra.{Double2Const, Expression, Zero} -import optimus.optimization.MPModel -import optimus.optimization.enums.{SolutionStatus, SolverLib} -import optimus.optimization.model.MPFloatVar -import org.slf4j.Logger -import squants.{Power, Time} - -import java.util.UUID - -/** Flex strategy that optimizes over a fixed amount of time steps into the - * future. Works with [[MathFlexOptions]], which describe constraints on - * present and future behavior of a participant. - * - * @param sampleTime - * The amount of time between the steps. - * @param predictionHorizon - * The amount of time that is predicted into the future, i.e. the last step - * is this amount of time away from now. - * @param powerObjectiveFactory - * The objective to optimize for. - * @param logger - * The logger to use. - */ -final case class OptimizedFlexStrat( - sampleTime: Time, - predictionHorizon: Time, - powerObjectiveFactory: PowerObjectiveFactory, - logger: Logger, -) extends EmModelStrat[MathFlexOptions[?, ? <: OperationVars]] { - - override def determineFlexControl( - flexOptions: Iterable[ - (? <: AssetInput, MathFlexOptions[?, ? <: OperationVars]) - ], - target: Power, - currentTick: Long, - ): Iterable[(UUID, Power)] = { - - given model: MPModel = MPModel(SolverLib.oJSolver) - - val sampleTicks = sampleTime.toSeconds.toLong - val lastPredictedTick = currentTick + predictionHorizon.toSeconds.toLong - - val ticks = Range.Long(currentTick, lastPredictedTick, sampleTicks) - - val assetVars = flexOptions.map { case (asset: AssetInput, fo) => - addAssetConstraints(asset.getUuid, fo, ticks) - } - - val objectiveContainer = - buildObjective(assetVars, target, sampleTime, powerObjectiveFactory) - - model.minimize(objectiveContainer.objective) - - model.start() - - if model.getStatus != SolutionStatus.OPTIMAL then - throw new CriticalFailureException( - s"Optimization ended with unexpected status ${model.getStatus}, ${SolutionStatus.OPTIMAL} was expected." - ) - - objectiveContainer.softConstraints - .filter(_.getError > softConstraintThreshold) - .foreach { constraint => - logger.warn(constraint.getWarningMessage) - } - - // we're only interested in the solutions for the current time step - val assetCtrl = assetVars.map { - case AssetVarContainer(assetUuid, _, operationVars) => - val firstOp = operationVars(0) - assetUuid -> firstOp.getPowerSolution.getOrElse( - throw new CriticalFailureException( - s"No solution present for operation variables ${firstOp.getPowerExpression}" - ) - ) - } - - model.release() - - assetCtrl - } - - override def adaptFlexOptions( - assetInput: AssetInput, - flexOptions: MathFlexOptions[?, ? <: OperationVars], - ): MathFlexOptions[?, ? <: OperationVars] = flexOptions -} - -object OptimizedFlexStrat { - - /** Threshold for the error of soft constraints after optimization. Every soft - * constraint error should stay below this threshold. - */ - private val softConstraintThreshold: Double = 1e-3 - - /** Creates and adds constraints for the given asset for all given sample - * times. States and operating points are strung together according to the - * sample times. - * - * @param assetUuid - * The UUID of the asset. - * @param flexOptions - * The flex options that the asset provided. - * @param ticks - * The ticks of the sample times to add constraints for. - * @param model - * The optimization model to use. - * @tparam SV - * The type of state variables. - * @tparam OV - * The type of operation variables. - * @return - * A container that holds all state and operation variables. - */ - def addAssetConstraints[SV, OV <: OperationVars]( - assetUuid: UUID, - flexOptions: MathFlexOptions[SV, OV], - ticks: Seq[Long], - )(using model: MPModel): AssetVarContainer[SV, OV] = { - - val firstTick = ticks.headOption.getOrElse( - throw new CriticalFailureException( - "No ticks to add constraints for were given." - ) - ) - val otherTicks = ticks.tail - - val initialState = flexOptions.addInitialState(firstTick) - - val (allStates, allOperationVars) = - otherTicks.foldLeft(IndexedSeq(initialState), IndexedSeq.empty[OV]) { - case ((states, operationVars), tick) => - val addOp = flexOptions.addOperationConstraints(states.last) - val addState = - flexOptions.addNewStateConstraints( - states.last, - addOp, - tick, - ) - - (states.appended(addState), operationVars.appended(addOp)) - } - - AssetVarContainer(assetUuid, allStates, allOperationVars) - } - - /** Builds an objective to minimize given the asset variables and an objective - * factory. - * - * @param assetVars - * The asset vars to optimize for. - * @param target - * The target power for each time step. - * @param sampleTime - * The sample time. - * @param powerObjectiveFactory - * The factor for the objective to optimize at every time step. - * @param model - * The optimization model to use. - * @return - * An [[ObjectiveContainer]] holding the objective and soft constraints. - */ - def buildObjective( - assetVars: Iterable[AssetVarContainer[?, ? <: OperationVars]], - target: Power, - sampleTime: Time, - powerObjectiveFactory: PowerObjectiveFactory, - )(using model: MPModel): ObjectiveContainer = { - // asset vars should all have the same amount of operation vars, - // since they should have all been created with the same sample time steps - val timeSteps = assetVars.headOption.map(_.operationVars.size).getOrElse(0) - - val (objectiveResult, softConstraintsResult) = - Range(0, timeSteps) - .map { timeStep => - assetVars.map { - _.operationVars(timeStep) - } - } - .foldLeft[(Expression, Seq[SoftConstraint])](Zero, Seq.empty) { - case ((objective, allConstraints), opVars) => - val difference = opVars.foldLeft[Expression](Zero) { - case (powers, op) => - powers + op.getPowerExpression - } - target.toKilowatts - - val constraints = - opVars.flatMap(_.getSoftConstraint(sampleTime)) - val constraintsExpression = constraints - .map(_.getExpression) - .reduceLeftOption(_ + _) - .getOrElse(Zero) - - val powerObjective = powerObjectiveFactory.build(difference) - - ( - objective + constraintsExpression + powerObjective, - allConstraints.appendedAll(constraints), - ) - } - - ObjectiveContainer(objectiveResult, softConstraintsResult) - } - - /** Trait for factories of power objectives. An objective is created for the - * sum of power of a single time step. - */ - trait PowerObjectiveFactory { - - /** Creates an objective for a single time step involving the sum of power. - * - * @param totalPower - * The sum of power of all assets for a time step. - * @param model - * The optimization model to use. - * @return - * The objective as an expression. - */ - def build(totalPower: Expression)(using model: MPModel): Expression - } - - /** Creates an objective that simply minimizes the absolute value of the sum - * of power by using an epigraph constraint. - */ - object MinAbsPowerObjectiveFactory extends PowerObjectiveFactory { - - override def build( - totalPower: Expression - )(using model: MPModel): Expression = { - val d = MPFloatVar.positive("d") - model.add(d >:= totalPower) - model.add(d >:= -totalPower) - - d - } - - } - - /** Creates an objective that uses a piecewise-linear (over-)approximation of - * the quadratic function on the sum of power. The convex epigraph is used to - * derive a linear constraint. Effectively, higher power values are punished - * more than lower ones. - * - * The piecewise approximation is created with a fixed number of segments - * (secant lines) up until given last segment. - * - * @param segmentCount - * The number of segments (secant lines) to create. A high number of - * segments might impact efficiency. - * @param lastSegment - * The value of the last segment boundary. This should be set close to the - * maximum value that is to be expected, otherwise the approximation - * becomes inaccurate. - */ - class LinearizedQuadraticPowerObjectiveFactory( - segmentCount: Int, - lastSegment: Double, - ) extends PowerObjectiveFactory { - - override def build( - totalPower: Expression - )(using model: MPModel): Expression = { - - val powerAbs = MPFloatVar.positive("powerAbs") - model.add(powerAbs >:= totalPower) - model.add(powerAbs >:= -totalPower) - - val segmentSize = lastSegment / segmentCount - - val t = MPFloatVar.positive("t") - - Range.inclusive(0, segmentCount).map(_ * segmentSize).sliding(2).foreach { - case Seq(uCurrent, uNext) => - val m = uCurrent + uNext - val b = -uCurrent * uNext - - model.add(t >:= m * powerAbs + b) - } - - // normalize the final value so that it maximizes - // somewhat close to the absolute value - val normalizationFactor = 1 / lastSegment - - t * normalizationFactor - } - - } - - /** Container holding all variables for one asset. - * - * @param assetUuid - * The UUID of the asset. - * @param states - * All state variables of the asset. - * @param operationVars - * The operation variables of the asset. - * @tparam SV - * The type of state variables. - * @tparam OV - * The type of operation variables. - */ - final case class AssetVarContainer[SV, OV <: OperationVars]( - assetUuid: UUID, - states: IndexedSeq[SV], - operationVars: IndexedSeq[OV], - ) - - /** Container holding the complete objective to minimize and relevant soft - * constraints. - * - * @param objective - * The objective, including all soft constraint expressions. - * @param softConstraints - * All soft constraints. - */ - final case class ObjectiveContainer( - objective: Expression, - softConstraints: Iterable[SoftConstraint], - ) - -} diff --git a/src/main/scala/edu/ie3/simona/model/em/opt/OptimizedFlexStrat.scala b/src/main/scala/edu/ie3/simona/model/em/opt/OptimizedFlexStrat.scala new file mode 100644 index 0000000000..c583188532 --- /dev/null +++ b/src/main/scala/edu/ie3/simona/model/em/opt/OptimizedFlexStrat.scala @@ -0,0 +1,445 @@ +/* + * © 2025. TU Dortmund University, + * Institute of Energy Systems, Energy Efficiency and Energy Economics, + * Research group Distribution grid planning and operation + */ + +package edu.ie3.simona.model.em.opt + +import edu.ie3.datamodel.models.input.AssetInput +import edu.ie3.simona.exceptions.CriticalFailureException +import edu.ie3.simona.model.em.EmModelStrat +import edu.ie3.simona.model.em.opt.OptimizedFlexStrat.* +import edu.ie3.simona.model.em.opt.SoftConstraint.AbsValueSoftConstraint +import edu.ie3.simona.ontology.messages.flex.EnergyBoundariesFlexOptions +import edu.ie3.simona.ontology.messages.flex.EnergyBoundariesFlexOptions.AssetEnergyBoundaries +import edu.ie3.util.interval.ClosedInterval +import edu.ie3.util.scala.quantities.DefaultQuantities.zeroKWh +import optimus.algebra.{Const, Double2Const, Expression, Zero} +import optimus.optimization.MPModel +import optimus.optimization.enums.{SolutionStatus, SolverLib} +import optimus.optimization.model.{MPFloatVar, MPVar} +import org.slf4j.Logger +import squants.energy.Kilowatts +import squants.{Each, Power, Time} + +import java.util.UUID + +/** Energy management strategy that optimizes flexibility usage over a fixed + * amount of time steps into the future. Takes [[EnergyBoundariesFlexOptions]] + * as inputs, which provide the required energy boundaries and power limits for + * constraining asset operation. + * + * @param sampleTime + * The amount of time between the steps. + * @param predictionHorizon + * The amount of time that is predicted into the future, i.e. the last step + * is this amount of time away from the current point in simulation time. + * Should be a multiple of [[sampleTime]]. + * @param powerObjectiveFactory + * A factory creating the optimization objective to use. + * @param logger + * The logger to use. + */ +final case class OptimizedFlexStrat( + sampleTime: Time, + predictionHorizon: Time, + powerObjectiveFactory: PowerObjectiveFactory, + logger: Logger, +) extends EmModelStrat[EnergyBoundariesFlexOptions] { + + override def determineFlexControl( + flexOptions: Iterable[(? <: AssetInput, EnergyBoundariesFlexOptions)], + target: Power, + currentTick: Long, + ): Iterable[(UUID, Power)] = { + + given model: MPModel = MPModel(SolverLib.oJSolver) + + val sampleTicks = sampleTime.toSeconds.toLong + val lastPredictedTick = currentTick + predictionHorizon.toSeconds.toLong + + val ticks = Range.Long(currentTick, lastPredictedTick, sampleTicks) + + val assetVars = addAssetConstraints( + flexOptions.map { case (asset: AssetInput, fo) => asset.getUuid -> fo }, + sampleTime, + ticks, + ) + + val objectiveContainer = + buildObjective(assetVars, target, powerObjectiveFactory) + + model.minimize(objectiveContainer.objective) + + model.start() + + if model.getStatus != SolutionStatus.OPTIMAL then + throw new CriticalFailureException( + s"Optimization ended with unexpected status ${model.getStatus}, ${SolutionStatus.OPTIMAL} was expected." + ) + + objectiveContainer.softConstraints + .filter(_.getError > softConstraintThreshold) + .foreach { constraint => + logger.warn(constraint.getWarningMessage) + } + + // we're only interested in the solutions for the current time step + val assetCtrl = assetVars.map { + case AssetVarContainer(assetUuid, results) => + val setPoint = results + .map { + // Taking only the first result for set points + _.headOption + .getOrElse( + throw new CriticalFailureException( + s"Empty results for asset $assetUuid" + ) + ) + // Operating point of first result + .getOperationResult + } + .reduceOption(_ + _) + .getOrElse( + throw new CriticalFailureException( + s"No results present for asset $assetUuid" + ) + ) + assetUuid -> setPoint + } + + model.release() + + assetCtrl + } + + override def adaptFlexOptions( + assetInput: AssetInput, + flexOptions: EnergyBoundariesFlexOptions, + ): EnergyBoundariesFlexOptions = flexOptions + +} + +object OptimizedFlexStrat { + + /** Threshold for the error of soft constraints after optimization. Every soft + * constraint error should stay below this threshold. + */ + private val softConstraintThreshold: Double = 1e-3 + + /** Creates and adds constraints for the given flex options for given sample + * times. States and operating points are linked according to the time steps. + * + * @param flexOptions + * The flex options that connected assets provided. + * @param sampleTime + * The amount of time between the steps. + * @param ticks + * The ticks (including current tick to last predicted tick) of the time + * steps to add constraints for. Should all be sample time duration apart + * from each other. + * @param model + * The optimization model to add variables and constraints to. + * @return + * Containers that holds all results, including state and operation + * variables. + */ + def addAssetConstraints( + flexOptions: Iterable[(UUID, EnergyBoundariesFlexOptions)], + sampleTime: Time, + ticks: Seq[Long], + )(using model: MPModel): Iterable[AssetVarContainer] = { + flexOptions.map { case (assetUUID, fo) => + val results = fo.energyBoundaries + .map(adaptEnergyBoundaries) + .map { boundaries => + ticks.tail.foldLeft[IndexedSeq[StepResults]](IndexedSeq.empty) { + case (previousResults, tick) => + val previousState = previousResults.lastOption.flatMap(_.state) + + val res = addAssetStep( + boundaries, + tick, + sampleTime, + previousState, + ) + + previousResults.appended(res) + } + } + + AssetVarContainer(assetUUID, results) + } + } + + /** Creates and adds constraints for a single asset and time step. Relevant + * variables and soft constraints are returned in a [[StepResults]] + * container. + * + * @param energyBoundaries + * The asset energy boundaries to use. + * @param tick + * The tick to add constraints for. + * @param sampleTime + * The amount of time between the steps. + * @param maybePreviousState + * Optionally, the previous state variable. + * @param model + * The optimization model to add variables and constraints to. + * @return + * The results for this asset and time step. + */ + private def addAssetStep( + energyBoundaries: AssetEnergyBoundaries, + tick: Long, + sampleTime: Time, + maybePreviousState: Option[Expression], + )(using model: MPModel): StepResults = { + + val energyLimits = energyBoundaries.energyLimits + .maxBefore(tick + 1) + .map { case (_, limits) => + limits + } + .getOrElse(throw new CriticalFailureException("No energy limits found!")) + + val formerEnergyLimits = + energyBoundaries.energyLimits.maxBefore(tick).map { case (_, limits) => + limits + } + + if energyLimits.getLower == energyLimits.getUpper && + formerEnergyLimits.forall(limits => limits.getLower == limits.getUpper) + then { + // there is no flexibility at all, thus we don't need any state to keep track of + + val formerEnergy = formerEnergyLimits.map(_.getUpper).getOrElse(zeroKWh) + val currentEnergy = energyLimits.getUpper + + val fixedPower = (currentEnergy - formerEnergy) / sampleTime + StepResults(Const(fixedPower.toKilowatts), None, None) + } else { + // we do have some flexibility at this point in time, model it + + // we use charging efficiency for both charging and discharging, + // since we use the adapted storage model here + val eta = energyBoundaries.etaCharge + + // determining a previous state + val previousState = maybePreviousState.getOrElse { + + // we have been given no former state as a parameter. Either... + formerEnergyLimits + // ... there was no flexibility in the last step, thus we use the last energy value + .filter(limits => limits.getLower == limits.getUpper) + .map(limits => Const(limits.getUpper.toKilowattHours)) + // ... or this is the initial step, thus we start at 0 + .getOrElse(Const(0d)) + } + + // modeling the operating point (power), + // valid between that previous and new state + val p = MPFloatVar( + symbol = "p", + lowerBound = energyBoundaries.powerLimits.getLower.toKilowatts, + upperBound = energyBoundaries.powerLimits.getUpper.toKilowatts, + ) + + // modeling the new state (stored energy) + val newState = + if energyLimits.getLower == energyLimits.getUpper then + Const(energyLimits.getUpper.toKilowattHours) + else + MPFloatVar( + symbol = "state", + lowerBound = energyLimits.getLower.toKilowattHours, + upperBound = energyLimits.getUpper.toKilowattHours, + ) + + val softConstraint = + if eta == Each(1) then { + // there are no charging/discharging losses, we can keep it simple + + model.add(newState := previousState + p * sampleTime.toHours) + None + } else { + // there are charging/discharging losses, thus use the full model + + val pAbsMax = + energyBoundaries.powerLimits.getUpper.max( + -energyBoundaries.powerLimits.getLower + ) + + val pAbs = MPFloatVar( + symbol = "pAbs", + lowerBound = 0, + upperBound = pAbsMax.toKilowatts, + ) + + model.add(pAbs >:= p) + model.add(pAbs >:= -p) + + model.add( + newState := previousState + (p - pAbs * (1 - eta.toEach)) * sampleTime.toHours + ) + + Some(AbsValueSoftConstraint(p, pAbs, eta, sampleTime)) + } + + StepResults(p, Some(newState), softConstraint) + } + + } + + /** Creates flex options that are equivalent to the original with regard to + * optimization, which are able to be optimized with a linear model though. + * In order to achieve this, a common efficiency needs to be calculated for + * charging and discharging operations, eliminating the need to distinguish + * between charging and discharging when formulating the state constraint. + * Furthermore, energy limits are adapted to work with the adapted + * efficiency. + * + * @param boundaries + * The energy boundaries to be adapted. + * @return + * The adapted energy boundaries. + */ + def adaptEnergyBoundaries( + boundaries: AssetEnergyBoundaries + ): AssetEnergyBoundaries = { + val etaCh = boundaries.etaCharge.toEach + val etaDis = boundaries.etaDischarge.toEach + + val etaAvg = (2 * etaCh * etaDis) / (1 + etaCh * etaDis) + + val newEnergyLimits = boundaries.energyLimits.map { case (tick, limits) => + val newLower = (limits.getLower / etaCh) * etaAvg + val newUpper = (limits.getUpper / etaCh) * etaAvg + + tick -> ClosedInterval(newLower, newUpper) + } + + val etaAvgEach = Each(etaAvg) + + boundaries.copy( + energyLimits = newEnergyLimits, + etaCharge = etaAvgEach, + etaDischarge = etaAvgEach, + ) + + } + + /** Builds an objective to minimize given the asset variables and an objective + * factory. + * + * @param assetVars + * The asset vars to optimize for. + * @param target + * The target power for each time step. + * @param powerObjectiveFactory + * The factor for the objective to optimize at every time step. + * @param model + * The optimization model to add variables and constraints to. + * @return + * An [[ObjectiveContainer]] holding the objective and soft constraints. + */ + def buildObjective( + assetVars: Iterable[AssetVarContainer], + target: Power, + powerObjectiveFactory: PowerObjectiveFactory, + )(using model: MPModel): ObjectiveContainer = { + // asset vars should all have the same amount of operation vars, + // since they should have all been created with the same sample time steps + val timeSteps = assetVars.headOption + .flatMap(_.results.headOption.map(_.size)) + .getOrElse(0) + + val (objectiveResult, softConstraintsResult) = + Range(0, timeSteps) + // first, sort all results by time step + .map { timeStep => + assetVars.flatMap { + _.results.map(_(timeStep)) + } + } + // then, build objective for every time step and combine them + .foldLeft[(Expression, Seq[SoftConstraint])](Zero, Seq.empty) { + case ((objective, allConstraints), results) => + val difference = results + .map(_.operation) + .reduceOption[Expression](_ + _) + .getOrElse(Zero) + val powerObjective = powerObjectiveFactory.build(difference, target) + + val constraints = + results.flatMap(_.softConstraint) + val constraintsExpression = constraints + .map(_.getExpression) + .reduceLeftOption(_ + _) + .getOrElse(Zero) + + ( + objective + constraintsExpression + powerObjective, + allConstraints.appendedAll(constraints), + ) + } + + ObjectiveContainer(objectiveResult, softConstraintsResult) + } + + /** Container holding all optimization steps (including variables and soft + * constraints) for one asset. + * + * @param assetUuid + * The UUID of the asset. + * @param results + * All step results of the asset, including state and operation variables. + */ + final case class AssetVarContainer( + assetUuid: UUID, + results: Seq[IndexedSeq[StepResults]], + ) + + /** Container holding the relevant variables and potentially a soft constraint + * for a specific asset and optimization time step. + * + * @param operation + * The operating point between the previous state and [[state]]. + * @param state + * Optionally the state that follows the operating point. + * @param softConstraint + * Optionally a soft constraint. + */ + final case class StepResults( + operation: Const | MPVar, + state: Option[Expression], + softConstraint: Option[SoftConstraint], + ) { + + def getOperationResult: Power = Kilowatts(operation match { + case const: Const => const.value + case variable: MPVar => + variable.value.getOrElse( + throw new CriticalFailureException( + s"No result present for variable $variable" + ) + ) + }) + + } + + /** Container holding the complete objective to minimize and relevant soft + * constraints. + * + * @param objective + * The objective, including all soft constraint expressions. + * @param softConstraints + * All soft constraints. + */ + final case class ObjectiveContainer( + objective: Expression, + softConstraints: Iterable[SoftConstraint], + ) + +} diff --git a/src/main/scala/edu/ie3/simona/model/em/opt/PowerObjectiveFactory.scala b/src/main/scala/edu/ie3/simona/model/em/opt/PowerObjectiveFactory.scala new file mode 100644 index 0000000000..23c1d7dfd5 --- /dev/null +++ b/src/main/scala/edu/ie3/simona/model/em/opt/PowerObjectiveFactory.scala @@ -0,0 +1,112 @@ +/* + * © 2025. TU Dortmund University, + * Institute of Energy Systems, Energy Efficiency and Energy Economics, + * Research group Distribution grid planning and operation + */ + +package edu.ie3.simona.model.em.opt + +import optimus.algebra.{Double2Const, Expression} +import optimus.optimization.MPModel +import optimus.optimization.model.{MPFloatVar, MPVar} +import squants.Power + +/** Trait for factories of power objectives. An objective is created for the sum + * of power of a single time step. + */ +trait PowerObjectiveFactory { + + /** Creates an objective expression involving the sum of power for a single + * time step. + * + * @param totalPower + * The sum of power of all assets for a time step. + * @param model + * The optimization model to add variables and constraints to. + * @return + * The objective as an expression. + */ + def build(totalPower: Expression, target: Power)(using + model: MPModel + ): Expression +} + +object PowerObjectiveFactory { + + /** Creates an objective that simply minimizes the absolute value of the sum + * of power by using an epigraph constraint. + */ + object MinAbsPowerObjectiveFactory extends PowerObjectiveFactory { + + override def build( + totalPower: Expression, + target: Power, + )(using model: MPModel): Expression = { + val difference = totalPower - target.toKilowatts + + absoluteValue(difference, "differenceAbs") + } + + } + + /** Creates an objective that uses a piecewise-linear (over-)approximation of + * the quadratic function on the sum of power. The convex epigraph is used to + * derive a linear constraint. Effectively, higher power values are punished + * more than lower ones. + * + * The piecewise approximation is created with a fixed number of segments + * (secant lines) up until given last segment. + * + * @param segmentCount + * The number of segments (secant lines) to create. Increasing the number + * of segments improves the accuracy of the approximation, but might impact + * efficiency. + * @param lastSegment + * The value of the last segment boundary. This should be set close to the + * maximum value that is to be expected, otherwise the approximation + * becomes inaccurate beyond this value. + */ + class LinearizedQuadraticPowerObjectiveFactory( + segmentCount: Int, + lastSegment: Double, + ) extends PowerObjectiveFactory { + + override def build( + totalPower: Expression, + target: Power, + )(using model: MPModel): Expression = { + val difference = totalPower - target.toKilowatts + + val differenceAbs = absoluteValue(difference, "differenceAbs") + + val segmentSize = lastSegment / segmentCount + + val t = MPFloatVar.positive("t") + + Range.inclusive(0, segmentCount).map(_ * segmentSize).sliding(2).foreach { + case Seq(uCurrent, uNext) => + val m = uCurrent + uNext + val b = -uCurrent * uNext + + model.add(t >:= m * differenceAbs + b) + } + + // normalize the final value so that it maximizes + // somewhat close to the absolute value + val normalizationFactor = 1 / lastSegment + + t * normalizationFactor + } + + } + + private def absoluteValue(value: Expression, name: String)(using + model: MPModel + ): MPVar = { + val valueAbs = MPFloatVar.positive(name) + model.add(valueAbs >:= value) + model.add(valueAbs >:= -value) + valueAbs + } + +} diff --git a/src/main/scala/edu/ie3/simona/model/em/opt/SoftConstraint.scala b/src/main/scala/edu/ie3/simona/model/em/opt/SoftConstraint.scala new file mode 100644 index 0000000000..0dd0a2e939 --- /dev/null +++ b/src/main/scala/edu/ie3/simona/model/em/opt/SoftConstraint.scala @@ -0,0 +1,125 @@ +/* + * © 2025. TU Dortmund University, + * Institute of Energy Systems, Energy Efficiency and Energy Economics, + * Research group Distribution grid planning and operation + */ + +package edu.ie3.simona.model.em.opt + +import edu.ie3.simona.exceptions.CriticalFailureException +import optimus.algebra.{Double2Const, Expression} +import optimus.optimization.model.MPFloatVar +import squants.{Dimensionless, Time} + +/** Trait to be extended by classes detailing a soft constraint as part of the + * optimization objective, including possible error handling. + */ +trait SoftConstraint { + + /** The soft constraint expression to be included in the objective to be + * minimized. + * + * @return + * The soft constraint expression. + */ + def getExpression: Expression + + /** Returns the amount of error stemming from the soft constraint. Only call + * this if you're sure that a solution has been determined! A + * [[edu.ie3.simona.exceptions.CriticalFailureException]] will be thrown + * otherwise. + * + * @return + * The amount of error. + */ + def getError: Double + + /** A warning message explaining what was expected and what happened instead. + * The message only makes sense if the error is actually larger than + * expected. + * + * @return + * The warning message. + */ + def getWarningMessage: String + +} + +object SoftConstraint { + + /** Soft constraint for an absolute value of a free variable. + * + * @param variable + * The variable that can be assigned a positive or negative number. + * @param absoluteVariable + * The variable that is supposed to be set to the absolute value of the + * [[variable]]. + * @param penalty + * The penalty factor to be multiplied with the absolute value for the soft + * constraint expression. + */ + final case class AbsValueSoftConstraint( + variable: MPFloatVar, + absoluteVariable: MPFloatVar, + penalty: Double, + ) extends SoftConstraint { + + override def getExpression: Expression = + absoluteVariable * penalty + + override def getError: Double = { + val (value, absoluteValue) = getVals + math.abs(math.abs(value) - absoluteValue) + } + + override def getWarningMessage: String = { + val (value, absoluteValue) = getVals + s"Soft constraint for storage: Approximated absolute value $absoluteValue " + + s"and actual absolute value ${math.abs(value)} are $getError apart." + } + + private def getVals: (Double, Double) = + variable.value + .zip(absoluteVariable.value) + .getOrElse( + throw new CriticalFailureException( + "Solution are expected to be determined at this point!" + ) + ) + + } + + object AbsValueSoftConstraint { + + /** Creates a soft constraint for the absolute value of a power value, used + * during energy optimization. + * + * @param p + * The power value variable that can be positive or negative. + * @param pAbs + * The variable that is supposed to hold the absolute value of the power. + * @param eta + * The charging/discharging efficiency. + * @param duration + * The duration for which the power is applied. + * @return + * The soft constraint. + */ + def apply( + p: MPFloatVar, + pAbs: MPFloatVar, + eta: Dimensionless, + duration: Time, + ): AbsValueSoftConstraint = { + // Total penalty is slightly larger than the model losses. + // Thus, the value of pAbs should be pushed down to the + // absolute of p. + val epsilon = 1e-6 + + val penalty = (1 - eta.toEach + epsilon) * duration.toHours + + AbsValueSoftConstraint(p, pAbs, penalty) + } + } + +} diff --git a/src/main/scala/edu/ie3/simona/model/participant/BmModel.scala b/src/main/scala/edu/ie3/simona/model/participant/BmModel.scala index 542e2b9ba9..4f075f08f2 100644 --- a/src/main/scala/edu/ie3/simona/model/participant/BmModel.scala +++ b/src/main/scala/edu/ie3/simona/model/participant/BmModel.scala @@ -20,6 +20,7 @@ import edu.ie3.simona.model.participant.ParticipantModel.{ ParticipantModelFactory, } import edu.ie3.simona.model.participant.control.QControl +import edu.ie3.simona.model.participant.flex.ParticipantInflexiblePowerLimitFlexModel import edu.ie3.simona.ontology.messages.flex.FlexType import edu.ie3.simona.service.Data.PrimaryData import edu.ie3.simona.service.Data.PrimaryData.ComplexPower diff --git a/src/main/scala/edu/ie3/simona/model/participant/FixedFeedInModel.scala b/src/main/scala/edu/ie3/simona/model/participant/FixedFeedInModel.scala index 29c53aab73..3ed94cc256 100644 --- a/src/main/scala/edu/ie3/simona/model/participant/FixedFeedInModel.scala +++ b/src/main/scala/edu/ie3/simona/model/participant/FixedFeedInModel.scala @@ -13,6 +13,7 @@ import edu.ie3.datamodel.models.result.system.{ } import edu.ie3.simona.model.participant.ParticipantModel.* import edu.ie3.simona.model.participant.control.QControl +import edu.ie3.simona.model.participant.flex.ParticipantInflexiblePowerLimitFlexModel import edu.ie3.simona.ontology.messages.flex.FlexType import edu.ie3.simona.service.Data.PrimaryData.{ ComplexPower, diff --git a/src/main/scala/edu/ie3/simona/model/participant/PowerSeriesMathFlexOptions.scala b/src/main/scala/edu/ie3/simona/model/participant/PowerSeriesMathFlexOptions.scala deleted file mode 100644 index 8ba9d4849c..0000000000 --- a/src/main/scala/edu/ie3/simona/model/participant/PowerSeriesMathFlexOptions.scala +++ /dev/null @@ -1,83 +0,0 @@ -/* - * © 2025. TU Dortmund University, - * Institute of Energy Systems, Energy Efficiency and Energy Economics, - * Research group Distribution grid planning and operation - */ - -package edu.ie3.simona.model.participant - -import edu.ie3.simona.exceptions.CriticalFailureException -import edu.ie3.simona.model.participant.PowerSeriesMathFlexOptions.* -import edu.ie3.simona.ontology.messages.flex.MathFlexOptions -import edu.ie3.simona.ontology.messages.flex.MathFlexOptions.{ - OperationVars, - SoftConstraint, -} -import optimus.algebra.{Const, Expression} -import optimus.optimization.MPModel -import squants.{Power, Time} -import squants.energy.Kilowatts - -import scala.collection.immutable.SortedMap - -/** Flex options for participants that follow a fixed trajectory of power - * values. - * - * @param powers - * The power values as a [[SortedMap]], thus powers in between keys can be - * extract with `maxBefore`. - */ -final case class PowerSeriesMathFlexOptions(powers: SortedMap[Long, Power]) - extends MathFlexOptions[PowerStateVars, PowerOperationVars] { - - override def addInitialState(tick: Long)(using - model: MPModel - ): PowerStateVars = PowerStateVars(tick) - - override def addOperationConstraints( - state: PowerStateVars - )(using model: MPModel): PowerOperationVars = { - val (_, power) = powers - .maxBefore(state.tick + 1) - .getOrElse( - throw new CriticalFailureException( - s"No power found for tick ${state.tick} in provided power set $powers" - ) - ) - PowerOperationVars(Const(power.toKilowatts)) - } - - override def addNewStateConstraints( - formerState: PowerStateVars, - op: PowerOperationVars, - tick: Long, - )(using model: MPModel): PowerStateVars = PowerStateVars(tick) - -} - -object PowerSeriesMathFlexOptions { - - /** Only stores the tick, no actual state variables required. - * - * @param tick - * The tick of the state. - */ - final case class PowerStateVars(tick: Long) - - /** Stores the power as a constant. - * - * @param power - * The power value in kW. - */ - final case class PowerOperationVars(power: Const) extends OperationVars { - - override def getPowerExpression: Expression = power - - override def getPowerSolution: Option[Power] = Some(Kilowatts(power.value)) - - override def getSoftConstraint(duration: Time): Option[SoftConstraint] = - None - - } - -} diff --git a/src/main/scala/edu/ie3/simona/model/participant/PvModel.scala b/src/main/scala/edu/ie3/simona/model/participant/PvModel.scala index 13d5e44e3c..b88d71db5c 100644 --- a/src/main/scala/edu/ie3/simona/model/participant/PvModel.scala +++ b/src/main/scala/edu/ie3/simona/model/participant/PvModel.scala @@ -19,10 +19,13 @@ import edu.ie3.simona.model.participant.ParticipantModel.{ OperationChangeIndicator, ParticipantModelFactory, } -import edu.ie3.simona.model.participant.PvModel.{PvState, RadiationData} +import edu.ie3.simona.model.participant.PvModel.* import edu.ie3.simona.model.participant.SolarIrradiationCalculation.* import edu.ie3.simona.model.participant.control.QControl -import edu.ie3.simona.model.participant.flex.PowerSeriesMathFlexModel +import edu.ie3.simona.model.participant.flex.{ + ParticipantInflexiblePowerLimitFlexModel, + ParticipantInflexibleEnergyLimitFlexModel, +} import edu.ie3.simona.ontology.messages.flex.FlexType import edu.ie3.simona.service.Data.PrimaryData.{ ComplexPower, @@ -82,7 +85,7 @@ class PvModel private ( override val flexModels: Map[FlexType, ParticipantFlexModel[PvState]] = Map( FlexType.PowerLimit -> ParticipantInflexiblePowerLimitFlexModel(this), - FlexType.MathProgramming -> PowerSeriesMathFlexModel( + FlexType.EnergyBoundaries -> ParticipantInflexibleEnergyLimitFlexModel( this, _.toStateSeries, ), diff --git a/src/main/scala/edu/ie3/simona/model/participant/WecModel.scala b/src/main/scala/edu/ie3/simona/model/participant/WecModel.scala index 64c78b2edd..3c85a4e92a 100644 --- a/src/main/scala/edu/ie3/simona/model/participant/WecModel.scala +++ b/src/main/scala/edu/ie3/simona/model/participant/WecModel.scala @@ -22,7 +22,10 @@ import edu.ie3.simona.model.participant.ParticipantModel.{ } import edu.ie3.simona.model.participant.WecModel.* import edu.ie3.simona.model.participant.control.QControl -import edu.ie3.simona.model.participant.flex.PowerSeriesMathFlexModel +import edu.ie3.simona.model.participant.flex.{ + ParticipantInflexiblePowerLimitFlexModel, + ParticipantInflexibleEnergyLimitFlexModel, +} import edu.ie3.simona.model.system.Characteristic import edu.ie3.simona.model.system.Characteristic.XYPair import edu.ie3.simona.ontology.messages.flex.FlexType @@ -72,7 +75,7 @@ class WecModel private ( override val flexModels: Map[FlexType, ParticipantFlexModel[WecState]] = Map( FlexType.PowerLimit -> ParticipantInflexiblePowerLimitFlexModel(this), - FlexType.MathProgramming -> PowerSeriesMathFlexModel( + FlexType.EnergyBoundaries -> ParticipantInflexibleEnergyLimitFlexModel( this, _.toStateSeries, ), diff --git a/src/main/scala/edu/ie3/simona/model/participant/flex/PowerSeriesMathFlexModel.scala b/src/main/scala/edu/ie3/simona/model/participant/flex/ParticipantInflexibleEnergyLimitFlexModel.scala similarity index 69% rename from src/main/scala/edu/ie3/simona/model/participant/flex/PowerSeriesMathFlexModel.scala rename to src/main/scala/edu/ie3/simona/model/participant/flex/ParticipantInflexibleEnergyLimitFlexModel.scala index 145b5af71b..f50fd3688b 100644 --- a/src/main/scala/edu/ie3/simona/model/participant/flex/PowerSeriesMathFlexModel.scala +++ b/src/main/scala/edu/ie3/simona/model/participant/flex/ParticipantInflexibleEnergyLimitFlexModel.scala @@ -6,21 +6,21 @@ package edu.ie3.simona.model.participant.flex -import edu.ie3.simona.model.participant.{ - ParticipantFlexModel, - ParticipantModel, - PowerSeriesMathFlexOptions, -} +import edu.ie3.simona.model.participant.{ParticipantFlexModel, ParticipantModel} import edu.ie3.simona.model.participant.ParticipantModel.{ ModelState, OperatingPoint, } -import edu.ie3.simona.ontology.messages.flex.FlexOptions +import edu.ie3.simona.ontology.messages.flex.EnergyBoundariesFlexOptions.AssetEnergyBoundaries +import edu.ie3.simona.ontology.messages.flex.{ + EnergyBoundariesFlexOptions, + FlexOptions, +} import scala.collection.immutable.SortedMap /** Flex model implementation for [[ParticipantModel]]s producing - * [[PowerSeriesMathFlexOptions]] based on a forecast of power values. + * [[EnergyBoundariesFlexOptions]] based on a forecast of power values. * * @param model * The participant model to create forecast series for. @@ -30,7 +30,7 @@ import scala.collection.immutable.SortedMap * @tparam S * The type of state of the participant model. */ -class PowerSeriesMathFlexModel[S <: ModelState]( +class ParticipantInflexibleEnergyLimitFlexModel[S <: ModelState]( model: ParticipantModel[?, S], determineStates: S => SortedMap[Long, S], ) extends ParticipantFlexModel[S] { @@ -42,7 +42,7 @@ class PowerSeriesMathFlexModel[S <: ModelState]( tick -> op.activePower } - PowerSeriesMathFlexOptions(powerMap) + EnergyBoundariesFlexOptions(AssetEnergyBoundaries(powerMap)) } } diff --git a/src/main/scala/edu/ie3/simona/model/participant/ParticipantInflexiblePowerLimitFlexModel.scala b/src/main/scala/edu/ie3/simona/model/participant/flex/ParticipantInflexiblePowerLimitFlexModel.scala similarity index 88% rename from src/main/scala/edu/ie3/simona/model/participant/ParticipantInflexiblePowerLimitFlexModel.scala rename to src/main/scala/edu/ie3/simona/model/participant/flex/ParticipantInflexiblePowerLimitFlexModel.scala index c19fa3c352..ac3edd86e1 100644 --- a/src/main/scala/edu/ie3/simona/model/participant/ParticipantInflexiblePowerLimitFlexModel.scala +++ b/src/main/scala/edu/ie3/simona/model/participant/flex/ParticipantInflexiblePowerLimitFlexModel.scala @@ -4,9 +4,10 @@ * Research group Distribution grid planning and operation */ -package edu.ie3.simona.model.participant +package edu.ie3.simona.model.participant.flex import edu.ie3.simona.model.participant.ParticipantModel.ModelState +import edu.ie3.simona.model.participant.{ParticipantFlexModel, ParticipantModel} import edu.ie3.simona.ontology.messages.flex.{ FlexOptions, PowerLimitFlexOptions, diff --git a/src/main/scala/edu/ie3/simona/model/participant/load/LoadModel.scala b/src/main/scala/edu/ie3/simona/model/participant/load/LoadModel.scala index 255da94e53..192d670d29 100644 --- a/src/main/scala/edu/ie3/simona/model/participant/load/LoadModel.scala +++ b/src/main/scala/edu/ie3/simona/model/participant/load/LoadModel.scala @@ -12,17 +12,14 @@ import edu.ie3.datamodel.models.result.system.{ SystemParticipantResult, } import edu.ie3.simona.config.RuntimeConfig.LoadRuntimeConfig -import edu.ie3.simona.model.participant.{ - ParticipantFlexModel, - ParticipantModel, - ParticipantInflexiblePowerLimitFlexModel, -} +import edu.ie3.simona.model.participant.{ParticipantFlexModel, ParticipantModel} import edu.ie3.simona.model.participant.ParticipantModel.{ ActivePowerOperatingPoint, ModelState, OperationChangeIndicator, ParticipantModelFactory, } +import edu.ie3.simona.model.participant.flex.ParticipantInflexiblePowerLimitFlexModel import edu.ie3.simona.ontology.messages.flex.FlexType import edu.ie3.simona.service.Data.PrimaryData.{ ComplexPower, diff --git a/src/main/scala/edu/ie3/simona/model/participant/storage/StorageEnergyBoundariesFlexModel.scala b/src/main/scala/edu/ie3/simona/model/participant/storage/StorageEnergyBoundariesFlexModel.scala new file mode 100644 index 0000000000..033921e3b7 --- /dev/null +++ b/src/main/scala/edu/ie3/simona/model/participant/storage/StorageEnergyBoundariesFlexModel.scala @@ -0,0 +1,34 @@ +/* + * © 2025. TU Dortmund University, + * Institute of Energy Systems, Energy Efficiency and Energy Economics, + * Research group Distribution grid planning and operation + */ + +package edu.ie3.simona.model.participant.storage + +import edu.ie3.simona.model.participant.ParticipantFlexModel +import edu.ie3.simona.model.participant.storage.StorageModel.StorageState +import edu.ie3.simona.ontology.messages.flex.EnergyBoundariesFlexOptions.AssetEnergyBoundaries +import edu.ie3.simona.ontology.messages.flex.{ + EnergyBoundariesFlexOptions, + FlexOptions, +} + +class StorageEnergyBoundariesFlexModel(private val model: StorageModel) + extends ParticipantFlexModel[ + StorageState + ] { + + override def determineFlexOptions(state: StorageState): FlexOptions = + EnergyBoundariesFlexOptions( + AssetEnergyBoundaries( + eStorage = model.eStorage, + currentEnergy = state.storedEnergy, + pMax = model.pMax, + etaCharge = model.eta, + etaDischarge = model.eta, + currentTick = state.tick, + ) + ) + +} diff --git a/src/main/scala/edu/ie3/simona/model/participant/storage/StorageMathFlexModel.scala b/src/main/scala/edu/ie3/simona/model/participant/storage/StorageMathFlexModel.scala deleted file mode 100644 index 3e91bb811d..0000000000 --- a/src/main/scala/edu/ie3/simona/model/participant/storage/StorageMathFlexModel.scala +++ /dev/null @@ -1,233 +0,0 @@ -/* - * © 2025. TU Dortmund University, - * Institute of Energy Systems, Energy Efficiency and Energy Economics, - * Research group Distribution grid planning and operation - */ - -package edu.ie3.simona.model.participant.storage - -import edu.ie3.simona.exceptions.CriticalFailureException -import edu.ie3.simona.model.participant.ParticipantFlexModel -import edu.ie3.simona.model.participant.storage.StorageMathFlexModel.StorageMathFlexOptions -import edu.ie3.simona.model.participant.storage.StorageModel.StorageState -import edu.ie3.simona.ontology.messages.flex.MathFlexOptions.{ - OperationVars, - SoftConstraint, -} -import edu.ie3.simona.ontology.messages.flex.{FlexOptions, MathFlexOptions} -import optimus.algebra.{Const, Double2Const, Expression, Long2Const} -import optimus.optimization.MPModel -import optimus.optimization.model.{MPFloatVar, MPVar} -import squants.{Dimensionless, Each, Time} -import squants.energy.{Energy, Kilowatts, Power} - -/** Flex model implementation for [[StorageModel]] producing - * [[MathFlexOptions]]. - * - * @param model - * The [[StorageModel]] to use parameters from. - */ -class StorageMathFlexModel(private val model: StorageModel) - extends ParticipantFlexModel[ - StorageState - ] { - - override def determineFlexOptions(state: StorageState): FlexOptions = - StorageMathFlexOptions.createAdaptedFlexOptions( - state.storedEnergy, - model.eStorage, - model.pMax, - model.eta, - model.eta, - ) -} - -object StorageMathFlexModel { - - /** Flex options for a storage that uses adapted energy and efficiency values - * in order to stay in linear programming. This model can be derived from any - * storage model using separate charging and discharging efficiencies (refer - * to [[StorageMathFlexOptions.createAdaptedFlexOptions]] for more - * information). Here, a common loss is subtracted from stored energy with - * every charging and discharging operation. - * - * For this model, the absolute value of the charging/discharging power is - * required. To achieve this in linear programming, an epigraph constraint on - * the power variable is used, which requires a soft constraint as part of - * the objective. - * - * @param currentEnergy - * The stored energy currently stored within [[StorageModel]], adapted for - * linearity. - * @param eStorage - * The storage capacity of the [[StorageModel]], adapted for linearity. - * @param pMax - * The maximum charging and discharging power. - * @param eta - * The efficiency of charging and discharging operations, adapted for - * linearity. - */ - class StorageMathFlexOptions( - val currentEnergy: Energy, - val eStorage: Energy, - val pMax: Power, - val eta: Dimensionless, - ) extends MathFlexOptions[StorageStateVars, StorageOperationVars] { - - override def addInitialState( - tick: Long - )(using model: MPModel): StorageStateVars = { - val storedEnergy = Const(currentEnergy.toKilowattHours) - StorageStateVars(storedEnergy, tick) - } - - override def addOperationConstraints(state: StorageStateVars)(using - model: MPModel - ): StorageOperationVars = { - val p = MPFloatVar("p", -pMax.toKilowatts, pMax.toKilowatts) - val pAbs = MPFloatVar("pAbs", 0, pMax.toKilowatts) - - model.add(pAbs >:= p) - model.add(pAbs >:= -p) - - StorageOperationVars(p, pAbs, eta) - } - - override def addNewStateConstraints( - formerState: StorageStateVars, - op: StorageOperationVars, - tick: Long, - )(using model: MPModel): StorageStateVars = { - - val storedEnergy = MPFloatVar("storedEnergy", 0, eStorage.toKilowattHours) - val timeInHours = (tick - formerState.tick) / 3600 - - model.add( - storedEnergy := formerState.storedEnergy + (op.p - op.pAbs * (1 - eta.toEach)) * timeInHours - ) - - StorageStateVars(storedEnergy, tick) - } - } - - object StorageMathFlexOptions { - - /** Creates an equivalent linear battery model using parameters of a regular - * model. A common efficiency needs to be calculated for charging and - * discharging operations. Furthermore, energy amounts need to be adapted. - * - * @param currentEnergy - * The currently stored energy. - * @param eStorage - * The storage capacity. - * @param pMax - * The maximum charging and discharging power. - * @param etaCharging - * The charging efficiency. - * @param etaDischarging - * The discharging efficiency. - * @return - * An adapted model. - */ - def createAdaptedFlexOptions( - currentEnergy: Energy, - eStorage: Energy, - pMax: Power, - etaCharging: Dimensionless, - etaDischarging: Dimensionless, - ): StorageMathFlexOptions = { - - val etaCh = etaCharging.toEach - val etaDis = etaDischarging.toEach - - val etaAvg = (2 * etaCh * etaDis) / (1 + etaCh * etaDis) - - val adaptedCurrentEnergy = (currentEnergy / etaCh) * etaAvg - val adaptedEStorage = (eStorage / etaCh) * etaAvg - - new StorageMathFlexOptions( - adaptedCurrentEnergy, - adaptedEStorage, - pMax, - Each(etaAvg), - ) - } - - } - - /** Relevant data related to a storage state. - * - * @param storedEnergy - * The amount of stored energy either as a constant ([[Const]]) or as a - * variable ([[MPVar]]). - * @param tick - * The tick of this state. - */ - final case class StorageStateVars(storedEnergy: Expression, tick: Long) - - /** Relevant data related to a storage operating point. - * - * @param p - * The charging/discharging power in kW. Positive means charging. - * @param pAbs - * The absolute power in kW. This should approximate [[p]] closely. - * @param eta - * The common charging/discharging efficiency. - */ - final case class StorageOperationVars( - p: MPVar, - pAbs: MPVar, - eta: Dimensionless, - ) extends OperationVars { - - override def getPowerExpression: Expression = - p - - override def getPowerSolution: Option[Power] = - p.value.map(Kilowatts.apply) - - override def getSoftConstraint(duration: Time): Option[SoftConstraint] = { - // putting a penalty on pAbs, so that it comes - // as close as possible to the absolute power - Some(new StorageSoftConstraint(duration)) - } - - /** Soft constraint for storage operation variables. Required for [[pAbs]] - * to closely approximate the absolute value of [[p]]. - * - * @param duration - * The sample time, thus charging/discharging duration. - */ - private class StorageSoftConstraint(duration: Time) extends SoftConstraint { - - override def getExpression: Expression = { - // Total penalty is slightly larger than the losses - // calculated by StorageMathFlexOptions. Thus, the - // value of pAbs should be pushed down to the absolute - // of p. - val epsilon = 1e-6 - pAbs * (1 - eta.toEach + epsilon) * duration.toHours - } - - override def getError: Double = { - val (pValue, pAbsValue) = getVals - math.abs(math.abs(pValue) - pAbsValue) - } - - override def getWarningMessage: String = { - val (pValue, pAbsValue) = getVals - s"Soft constraint for storage: Approximated absolute power value $pAbsValue and absolute power value |$pValue| are $getError apart." - } - - private def getVals: (Double, Double) = p.value - .zip(pAbs.value) - .getOrElse( - throw new CriticalFailureException( - "Solution are expected to be determined at this point!" - ) - ) - - } - } - -} diff --git a/src/main/scala/edu/ie3/simona/model/participant/storage/StorageModel.scala b/src/main/scala/edu/ie3/simona/model/participant/storage/StorageModel.scala index 92979f831b..b69f7c6bc1 100644 --- a/src/main/scala/edu/ie3/simona/model/participant/storage/StorageModel.scala +++ b/src/main/scala/edu/ie3/simona/model/participant/storage/StorageModel.scala @@ -115,7 +115,7 @@ class StorageModel private ( override val flexModels: Map[FlexType, ParticipantFlexModel[StorageState]] = Map( FlexType.PowerLimit -> StoragePowerLimitFlexModel(this), - FlexType.MathProgramming -> StorageMathFlexModel(this), + FlexType.EnergyBoundaries -> StorageEnergyBoundariesFlexModel(this), ) override def determineState( diff --git a/src/main/scala/edu/ie3/simona/ontology/messages/flex/EnergyBoundariesFlexOptions.scala b/src/main/scala/edu/ie3/simona/ontology/messages/flex/EnergyBoundariesFlexOptions.scala new file mode 100644 index 0000000000..5590e290fc --- /dev/null +++ b/src/main/scala/edu/ie3/simona/ontology/messages/flex/EnergyBoundariesFlexOptions.scala @@ -0,0 +1,179 @@ +/* + * © 2025. TU Dortmund University, + * Institute of Energy Systems, Energy Efficiency and Energy Economics, + * Research group Distribution grid planning and operation + */ + +package edu.ie3.simona.ontology.messages.flex + +import edu.ie3.simona.exceptions.CriticalFailureException +import edu.ie3.simona.ontology.messages.flex.EnergyBoundariesFlexOptions.AssetEnergyBoundaries +import edu.ie3.util.interval.ClosedInterval +import edu.ie3.util.scala.quantities.DefaultQuantities.{zeroKW, zeroKWh} +import squants.time.Seconds +import squants.{Dimensionless, Each, Energy, Power} + +import scala.collection.immutable.SortedMap + +/** Energy boundaries for one or several assets. See [[AssetEnergyBoundaries]] + * for more details. + * + * @param energyBoundaries + * The energy boundaries. + */ +final case class EnergyBoundariesFlexOptions( + energyBoundaries: Seq[AssetEnergyBoundaries] +) extends FlexOptions + +object EnergyBoundariesFlexOptions { + + /** Creates energy boundaries with a single [[AssetEnergyBoundaries]]. + * + * @param singleBoundaries + * The [[AssetEnergyBoundaries]]. + * @return + */ + def apply( + singleBoundaries: AssetEnergyBoundaries + ): EnergyBoundariesFlexOptions = + EnergyBoundariesFlexOptions(Seq(singleBoundaries)) + + /** Energy boundaries for an asset. The energy limits (valid for the interval + * from tick to the next) constitute the boundaries between which flexibility + * can be used. + * + * @param energyLimits + * Energy limits that signify the potential upwards and downwards + * flexibility potential for the respective tick. The energy limits for all + * ticks relate to the energy potential at the current tick (which is + * defined to be zero). + * @param powerLimits + * The power limits, which limit the power of the complete asset for all + * time steps. If energy limits (upper and lower) are the same at some time + * step, power limits are ignored. + * @param etaCharge + * The charging efficiency. + * @param etaDischarge + * The discharging efficiency. + * @param tickDisconnect + * Optionally, the tick at which the storage will be disconnected, thus the + * upward or downward energy potential can not be used beyond this tick. + */ + final case class AssetEnergyBoundaries( + energyLimits: SortedMap[Long, ClosedInterval[Energy]], + powerLimits: ClosedInterval[Power], + etaCharge: Dimensionless = Each(1), + etaDischarge: Dimensionless = Each(1), + tickDisconnect: Option[Long] = None, + ) + + object AssetEnergyBoundaries { + + /** Creating energy boundaries for a fixed power time series. Assumes + * equidistant power series entries - otherwise, results are not defined! + * + * @param powerSeries + * The power time series (at equidistant ticks). Has to consist of at + * least two entries. + * @return + * The [[AssetEnergyBoundaries]]. + */ + def apply( + powerSeries: SortedMap[Long, Power] + ): AssetEnergyBoundaries = { + + val (firstTick, firstPower) = powerSeries.headOption.getOrElse( + throw new CriticalFailureException("Empty power time series!") + ) + + // adding a dummy tick at the end, so that the last actual power entry + // in the power time series is used + val lastSeriesTick = powerSeries.lastOption + .map { case (tick, _) => tick } + .getOrElse( + throw new CriticalFailureException("Empty power time series!") + ) + // dummy tick is the next logical tick of the tick series + // (e.g. (5, 10, 15, 20) -> 25) + val dummyTick = + lastSeriesTick + (lastSeriesTick - firstTick) / (powerSeries.size - 1) + + val (energySeries, _) = + powerSeries + // excluding first data, which we already extracted above + .tail + // adding dummy tick so that last actual power entry is recognized + .updated(dummyTick, zeroKW) + .foldLeft( + (SortedMap(firstTick -> zeroKWh), firstPower) + ) { case ((previousSeries, previousPower), (tick, power)) => + val (lastTick, lastEnergy) = previousSeries + .maxBefore(tick) + .getOrElse( + throw new CriticalFailureException( + s"No value before $tick in previous values $previousSeries" + ) + ) + + // added energy from last to current tick + val addedEnergy = previousPower * Seconds(tick - lastTick) + // total current energy + val tickEnergy = lastEnergy + addedEnergy + + (previousSeries.updated(tick, tickEnergy), power) + } + + val minPower = powerSeries.values.minOption.getOrElse(zeroKW) + val maxPower = powerSeries.values.maxOption.getOrElse(zeroKW) + + AssetEnergyBoundaries( + energyLimits = energySeries.map { case (tick, energy) => + tick -> ClosedInterval(energy, energy) + }, + powerLimits = ClosedInterval(minPower, maxPower), + ) + } + + /** Creating energy boundaries for a storage model with symmetrical maximum + * power. + * + * @param eStorage + * The storage capacity. + * @param currentEnergy + * The currently stored energy. + * @param pMax + * Maximum permissible active power. + * @param etaCharge + * Efficiency of the charging process. + * @param etaDischarge + * Efficiency of the discharging process. + * @param currentTick + * The current tick. + * @return + * [[AssetEnergyBoundaries]] for the storage model. + */ + def apply( + eStorage: Energy, + currentEnergy: Energy, + pMax: Power, + etaCharge: Dimensionless, + etaDischarge: Dimensionless, + currentTick: Long, + ): AssetEnergyBoundaries = + AssetEnergyBoundaries( + energyLimits = SortedMap( + currentTick -> new ClosedInterval( + // the downward energy potential is the currently stored energy + -currentEnergy, + // the upward energy potential is the energy to be stored until maximum capacity + eStorage - currentEnergy, + ) + ), + powerLimits = new ClosedInterval(-pMax, pMax), + etaCharge = etaCharge, + etaDischarge = etaDischarge, + ) + + } + +} diff --git a/src/main/scala/edu/ie3/simona/ontology/messages/flex/FlexType.scala b/src/main/scala/edu/ie3/simona/ontology/messages/flex/FlexType.scala index 73a0b6f601..9b3cecc652 100644 --- a/src/main/scala/edu/ie3/simona/ontology/messages/flex/FlexType.scala +++ b/src/main/scala/edu/ie3/simona/ontology/messages/flex/FlexType.scala @@ -10,5 +10,5 @@ package edu.ie3.simona.ontology.messages.flex * in subclasses of [[FlexOptions]]. */ enum FlexType { - case PowerLimit, MathProgramming + case PowerLimit, EnergyBoundaries } diff --git a/src/main/scala/edu/ie3/simona/ontology/messages/flex/MathFlexOptions.scala b/src/main/scala/edu/ie3/simona/ontology/messages/flex/MathFlexOptions.scala deleted file mode 100644 index 7cb13c30b3..0000000000 --- a/src/main/scala/edu/ie3/simona/ontology/messages/flex/MathFlexOptions.scala +++ /dev/null @@ -1,131 +0,0 @@ -/* - * © 2025. TU Dortmund University, - * Institute of Energy Systems, Energy Efficiency and Energy Economics, - * Research group Distribution grid planning and operation - */ - -package edu.ie3.simona.ontology.messages.flex - -import edu.ie3.simona.ontology.messages.flex.MathFlexOptions.OperationVars -import optimus.algebra.Expression -import optimus.optimization.MPModel -import squants.{Power, Time} - -/** Flex options used for mathematical programming, e.g. for optimization. - * Methods for creating a chain of states and operating points need to be - * provided. - * - * @tparam SV - * The type of state variables. - * @tparam OV - * The type of operation variables. - */ -trait MathFlexOptions[SV, OV <: OperationVars] extends FlexOptions { - - /** Adds and returns the initial state. - * - * @param tick - * The current tick, which is also the tick of the initial state. - * @param model - * The model to use. - * @return - * The initial state variable. - */ - def addInitialState(tick: Long)(using model: MPModel): SV - - /** Adds and returns operation variables for given state. - * - * @param state - * The state to add operation variables for. - * @param model - * The model to use. - * @return - * The operation variables. - */ - def addOperationConstraints(state: SV)(using model: MPModel): OV - - /** Adds and returns state variables for the state that results from the - * former state and given operation constraints. - * - * @param formerState - * The former state. - * @param op - * The operation variables and constraints. - * @param tick - * The tick to create the new state for. - * @param model - * The model to use. - * @return - * The state variable. - */ - def addNewStateConstraints(formerState: SV, op: OV, tick: Long)(using - model: MPModel - ): SV - -} - -object MathFlexOptions { - - /** Trait that needs to be extended by all types of operation variables. - */ - trait OperationVars { - - /** The final power expression to use within the objective of optimization. - * - * @return - * The power expression. - */ - def getPowerExpression: Expression - - /** The solution found for the power expression. Only available once - * optimization has run and succeeded. Will return [[None]] otherwise. - * - * @return - * The solution to the power expression, if applicable. - */ - def getPowerSolution: Option[Power] - - /** Returns the soft constraint for the operation variables, if applicable. - * - * @param duration - * The duration that the system participant was operating at this - * operating point, i.e. the sample time. - */ - def getSoftConstraint(duration: Time): Option[SoftConstraint] - - } - - /** Trait to be extended by classes detailing a soft constraint as part of the - * optimization objective and possible error handling. - */ - trait SoftConstraint { - - /** The soft constraint expression to be included in the objective to be - * minimized. - * - * @return - * The soft constraint expression. - */ - def getExpression: Expression - - /** Returns the amount of error stemming from the soft constraint. Only call - * this if you're sure that a solution has been determined! A - * [[edu.ie3.simona.exceptions.CriticalFailureException]] will be thrown - * otherwise. - * - * @return - * The amount of error. - */ - def getError: Double - - /** A warning message explaining what was expected and what happened - * instead. Only makes sense if the error is larger than expected. - * - * @return - * The warning message. - */ - def getWarningMessage: String - - } - -} diff --git a/src/test/scala/edu/ie3/simona/model/em/OptimizedFlexStratIT.scala b/src/test/scala/edu/ie3/simona/model/em/OptimizedFlexStratIT.scala deleted file mode 100644 index a5b97c58cd..0000000000 --- a/src/test/scala/edu/ie3/simona/model/em/OptimizedFlexStratIT.scala +++ /dev/null @@ -1,266 +0,0 @@ -/* - * © 2025. TU Dortmund University, - * Institute of Energy Systems, Energy Efficiency and Energy Economics, - * Research group Distribution grid planning and operation - */ - -package edu.ie3.simona.model.em - -import edu.ie3.simona.model.em.OptimizedFlexStrat.{ - AssetVarContainer, - LinearizedQuadraticPowerObjectiveFactory, - MinAbsPowerObjectiveFactory, -} -import edu.ie3.simona.model.participant.PowerSeriesMathFlexOptions -import edu.ie3.simona.model.participant.PowerSeriesMathFlexOptions.{ - PowerOperationVars, - PowerStateVars, -} -import edu.ie3.simona.model.participant.storage.StorageMathFlexModel.{ - StorageMathFlexOptions, - StorageOperationVars, - StorageStateVars, -} -import edu.ie3.simona.test.common.{MathFlexTestLike, UnitSpec} -import optimus.optimization.MPModel -import optimus.optimization.enums.{SolutionStatus, SolverLib} -import squants.Each -import squants.energy.{KilowattHours, Kilowatts} -import squants.time.Hours - -import java.util.UUID - -class OptimizedFlexStratIT extends UnitSpec with MathFlexTestLike { - - // Testing tolerances - given Double = 1e-6 - val constraintTolerance = 1e-3 - - "An optimized flex strat" when { - "provided with battery and constant constraints" should { - - given ticks: Seq[Long] = Range.Long.inclusive(0, 12 * 3600, 3600) - - // 33 kWh of feed-in in total, more than battery can store - def createPvVars(using - MPModel - ): AssetVarContainer[PowerStateVars, PowerOperationVars] = - OptimizedFlexStrat.addAssetConstraints( - UUID.fromString("0-0-0-0-1"), - PowerSeriesMathFlexOptions( - Seq(0, -6, -8, -7, -12, 0, 0, 0, 0, 0, 0, 0).toPowerMap - ), - ticks, - ) - - // 31 kWh of load in total, more than battery can provide - def createLoadVars(using - MPModel - ): AssetVarContainer[PowerStateVars, PowerOperationVars] = - OptimizedFlexStrat.addAssetConstraints( - UUID.fromString("0-0-0-0-2"), - PowerSeriesMathFlexOptions( - Seq(0, 0, 0, 0, 0, 0, 5, 10, 3, 7, 6, 0).toPowerMap - ), - ticks, - ) - - val batFo = StorageMathFlexOptions.createAdaptedFlexOptions( - KilowattHours(0), - KilowattHours(20), - Kilowatts(10), - Each(0.8), - Each(0.8), - ) - def createBatVars(using - MPModel - ): AssetVarContainer[StorageStateVars, StorageOperationVars] = - OptimizedFlexStrat.addAssetConstraints( - UUID.fromString("0-0-0-0-3"), - batFo, - ticks, - ) - - // since energy values have been adapted, we need this - // factor to convert back to "real" values - given EnergyConversionFactor = - EnergyConversionFactor(Each(0.8), batFo.eta) - - "compensate the load and feed-in with battery flexibility" in { - - given model: MPModel = MPModel(SolverLib.oJSolver) - - val pvVars = createPvVars - val loadVars = createLoadVars - val batVars = createBatVars - - val objectiveContainer = OptimizedFlexStrat.buildObjective( - Iterable(pvVars, loadVars, batVars), - Kilowatts(0), - Hours(1), - MinAbsPowerObjectiveFactory, - ) - - model.minimize(objectiveContainer.objective) - model.start(timeLimit = 10000) - - model.getStatus shouldBe SolutionStatus.OPTIMAL - - /* - EXPECTED RESULTS - Since excess power costs the same at all points in time and - at all magnitudes, there are many optimal solutions. - Thus, we only test for things that are true for every optimal - solution: We know when the battery should be definitely - full/empty and how much energy was charged/discharged. - */ - - { - objectiveContainer.softConstraints.foreach { constraint => - withClue(constraint.getWarningMessage) { - constraint.getError should be < constraintTolerance - } - } - - batVars.states(0).energyVal should approximate(0d) - - batVars.states.slice(0, 5).foreach { - _.energyVal should (be >= 0d and be < 20d) - } - - batVars.states.slice(5, 7).foreach { - _.energyVal should approximate(20d) - } - - batVars.states.slice(7, 11).foreach { - _.energyVal should (be >= 0d and be < 20d) - } - - batVars.states.slice(11, 13).foreach { - _.energyVal should approximate(0d) - } - - // we should've charged 20 kWh plus 5 kWh losses - val inputCharged = batVars.operationVars.slice(0, 6).map(_.pVal).sum - inputCharged should approximate(25) - - // we should've discharged 20 kWh minus 4 kWh losses - val outputDischarged = - batVars.operationVars.slice(6, 12).map(_.pVal).sum - outputDischarged should approximate(-16d) - - } withClue buildDebugString(batVars) - - model.release() - - } - - "minimize peaks" in { - - given model: MPModel = MPModel(SolverLib.oJSolver) - - val pvVars = createPvVars - val loadVars = createLoadVars - val batVars = createBatVars - - val objectiveContainer = OptimizedFlexStrat.buildObjective( - Iterable(pvVars, loadVars, batVars), - Kilowatts(0), - Hours(1), - LinearizedQuadraticPowerObjectiveFactory( - segmentCount = 10, - lastSegment = 10d, // 10 kW - ), - ) - - model.minimize(objectiveContainer.objective) - model.start(timeLimit = 10000) - - model.getStatus shouldBe SolutionStatus.OPTIMAL - - /* - EXPECTED RESULTS - When using the (linearized) quadratic objective, the - test is designed in a way such that only one optimal - solution exists: During the first phase (charging), - battery power is utilized so that 2 kW remains at every - time step. During the second phase (discharging), there's - exactly enough energy available to go down to 3 kW at - every time step. - */ - - { - objectiveContainer.softConstraints.foreach { constraint => - withClue(constraint.getWarningMessage) { - constraint.getError should be < constraintTolerance - } - } - - batVars.states(0).energyVal should approximate(0d) - - // 0 kW to compensate - batVars.operationVars(0).pVal should approximate(0d) - batVars.states(1).energyVal should approximate(0d) - - // 6 kW of feed-in to compensate, 2 kW remains - batVars.operationVars(1).pVal should approximate(4d) - batVars.states(2).energyVal should approximate(3.2d) - - // 8 kW of feed-in to compensate, 2 kW remains - batVars.operationVars(2).pVal should approximate(6d) - batVars.states(3).energyVal should approximate(8d) - - // 7 kW of feed-in to compensate, 2 kW remains - batVars.operationVars(3).pVal should approximate(5d) - batVars.states(4).energyVal should approximate(12d) - - // 12 kW of feed-in to compensate, 2 kW remains - batVars.operationVars(4).pVal should approximate(10d) - batVars.states(5).energyVal should approximate(20d) - - // 0 kW to compensate - batVars.operationVars(5).pVal should approximate(0d) - batVars.states(6).energyVal should approximate(20d) - - // 5 kW of load to compensate, 3 kW remains - batVars.operationVars(6).pVal should approximate(-2d) - batVars.states(7).energyVal should approximate(17.5d) - - // 10 kW of load to compensate, 3 kW remains - batVars.operationVars(7).pVal should approximate(-7d) - batVars.states(8).energyVal should approximate(8.75d) - - // 3 kW of load to compensate, 3 kW remains - batVars.operationVars(8).pVal should approximate(0d) - batVars.states(9).energyVal should approximate(8.75d) - - // 7 kW of load to compensate, 3 kW remains - batVars.operationVars(9).pVal should approximate(-4d) - batVars.states(10).energyVal should approximate(3.75d) - - // 6 kW of load to compensate, 3 kW remains - batVars.operationVars(10).pVal should approximate(-3d) - batVars.states(11).energyVal should approximate(0d) - - // 0 kW to compensate - batVars.operationVars(11).pVal should approximate(0d) - batVars.states(12).energyVal should approximate(0d) - - // we should've charged 20 kWh plus 5 kWh losses - val inputCharged = batVars.operationVars.slice(0, 6).map(_.pVal).sum - inputCharged should approximate(25) - - // we should've discharged 20 kWh minus 4 kWh losses - val outputDischarged = - batVars.operationVars.slice(6, 12).map(_.pVal).sum - outputDischarged should approximate(-16d) - - } withClue buildDebugString(batVars) - - model.release() - - } - } - } - -} diff --git a/src/test/scala/edu/ie3/simona/model/em/opt/OptimizedFlexStratSpec.scala b/src/test/scala/edu/ie3/simona/model/em/opt/OptimizedFlexStratSpec.scala new file mode 100644 index 0000000000..00670ef09d --- /dev/null +++ b/src/test/scala/edu/ie3/simona/model/em/opt/OptimizedFlexStratSpec.scala @@ -0,0 +1,929 @@ +/* + * © 2025. TU Dortmund University, + * Institute of Energy Systems, Energy Efficiency and Energy Economics, + * Research group Distribution grid planning and operation + */ + +package edu.ie3.simona.model.em.opt + +import edu.ie3.simona.model.em.opt.OptimizedFlexStratSpec.* +import edu.ie3.simona.model.em.opt.PowerObjectiveFactory.{ + LinearizedQuadraticPowerObjectiveFactory, + MinAbsPowerObjectiveFactory, +} +import edu.ie3.simona.ontology.messages.flex.EnergyBoundariesFlexOptions +import edu.ie3.simona.ontology.messages.flex.EnergyBoundariesFlexOptions.AssetEnergyBoundaries +import edu.ie3.simona.test.common.{OptimizingTestLike, UnitSpec} +import edu.ie3.util.scala.quantities.DefaultQuantities.* +import optimus.optimization.MPModel +import optimus.optimization.enums.{SolutionStatus, SolverLib} +import org.scalatest.OptionValues +import squants.energy.{Energy, KilowattHours, Kilowatts, WattHours} +import squants.time.Hours +import squants.{Dimensionless, Each, Power, Time} + +import java.util.UUID + +class OptimizedFlexStratSpec extends UnitSpec with OptimizingTestLike { + + private val pvUUID = UUID.fromString("0-0-0-0-1") + private val loadUUID = UUID.fromString("0-0-0-0-2") + private val batUUID = UUID.fromString("0-0-0-0-3") + + private val sampleTime = Hours(1) + private val sampleTicks = sampleTime.toSeconds.toLong + + // Testing tolerances + given Double = 1e-6 + given Energy = WattHours(1e-9) + private val constraintTolerance = 1e-3 + + "An optimized flex strat" when { + + "provided with a flex energy model to adapt" should { + + "create an adapted model correctly" in { + + val currentEnergy = KilowattHours(10) + val eStorage = KilowattHours(20) + val pMax = Kilowatts(10) + + val etas = Seq(.6, .65, .7, .75, .8, .85, 0.9, .92, .95, .98, 1) + + forEvery(Table("etaCharging", etas*)) { etaCharging => + forEvery(Table("etaDischarging", etas*)) { etaDischarging => + + val classic = ClassicModel( + currentEnergy = currentEnergy, + eStorage = eStorage, + pMax = pMax, + etaCharging = Each(etaCharging), + etaDischarging = Each(etaDischarging), + ) + + // using the implementation to derive an adapted model + val adaptedFlexOptions = OptimizedFlexStrat.adaptEnergyBoundaries( + AssetEnergyBoundaries( + eStorage = eStorage, + currentEnergy = currentEnergy, + pMax = pMax, + etaCharge = Each(etaCharging), + etaDischarge = Each(etaDischarging), + 0L, + ) + ) + val adaptedLimits = adaptedFlexOptions.energyLimits(0) + + val adapted = AdaptedModel( + currentEnergy = -adaptedLimits.getLower, + eStorage = adaptedLimits.getUpper - adaptedLimits.getLower, + pMax = adaptedFlexOptions.powerLimits.getUpper, + etaAvg = adaptedFlexOptions.etaCharge, + ) + + // charging until full, with maximum power + val power1 = pMax + val duration1 = (eStorage - currentEnergy) / (power1 * etaCharging) + // discharging three quarters, with half power + val power2 = -pMax / 2 + val duration2 = eStorage * 0.75 / (-power2 * 1 / etaDischarging) + // charging until half, with quarter power + val power3 = pMax / 4 + val duration3 = eStorage / 4 / (power3 * etaCharging) + // discharging until empty, with maximum power + val power4 = -pMax + val duration4 = eStorage / 2 / (-power4 * 1 / etaDischarging) + + val results = IndexedSeq(classic, adapted).map( + _.charge(power1, duration1) + .charge(power2, duration2) + .charge(power3, duration3) + .charge(power4, duration4) + .currentEnergy + ) + + // comparing both models: + // battery should be exactly empty + results(0) should approximate(zeroKWh) + results(1) should approximate(zeroKWh) + + } + } + } + + } + + "provided with simple battery flex options with losses" should { + + given ticks: Seq[Long] = + Range.Long.inclusive(0, 4 * sampleTicks, sampleTicks) + + // low efficiency for simplicity of the test + val batFlex = EnergyBoundariesFlexOptions( + AssetEnergyBoundaries( + eStorage = KilowattHours(24), + currentEnergy = KilowattHours(12), + pMax = Kilowatts(10), + etaCharge = Each(0.8), + etaDischarge = Each(0.8), + currentTick = 0L, + ) + ) + + // since energy values have been adapted, we need this + // factor to convert back to "real" values + given EnergyConversionFactor = + EnergyConversionFactor( + batFlex.energyBoundaries.headOption.value.etaCharge, + OptimizedFlexStrat + .adaptEnergyBoundaries(batFlex.energyBoundaries.headOption.value) + .etaCharge, + ) + + "balance out additional power within maximum battery power" in { + + given model: MPModel = MPModel(SolverLib.oJSolver) + + val constFlex = EnergyBoundariesFlexOptions( + AssetEnergyBoundaries( + Seq(5, -10, 10, -2).toPowerMap + ) + ) + + val flexOptions = Map( + loadUUID -> constFlex, + batUUID -> batFlex, + ) + + val assetVars = + OptimizedFlexStrat.addAssetConstraints(flexOptions, sampleTime, ticks) + + assetVars.toSeq should have length 2 + assetVars.foreach(_.results should have length 1) + assetVars.foreach(_.results.foreach(_ should have length 4)) + + val objectiveContainer = OptimizedFlexStrat.buildObjective( + assetVars, + zeroKW, + MinAbsPowerObjectiveFactory, + ) + + model.minimize(objectiveContainer.objective) + model.start(timeLimit = 10000) + + model.getStatus shouldBe SolutionStatus.OPTIMAL + + /* + EXPECTED RESULTS + Battery should be able to fully cover the additional power + */ + + val batVars = assetVars + .find(_.assetUuid == batUUID) + .getOrElse(fail(s"No asset variables for battery ($batUUID) found.")) + + val batRes = batVars.results.headOption + .getOrElse(fail(s"Empty results for battery ($batUUID).")) + + { + objectiveContainer.softConstraints.foreach { constraint => + withClue(constraint.getWarningMessage) { + constraint.getError should be < constraintTolerance + } + } + + // discharging 5 kWh plus 1.25 kWh losses + batRes(0).pVal should approximate(-5) + batRes(0).energyVal should approximate(-6.25) + + // charging 10 kWh minus 2 kWh losses + batRes(1).pVal should approximate(10) + batRes(1).energyVal should approximate(1.75) + + // discharging 10 kWh plus 2.5 kWh losses + batRes(2).pVal should approximate(-10) + batRes(2).energyVal should approximate(-10.75) + + // charging 2 kWh minus 0.4 kWh losses + batRes(3).pVal should approximate(2) + batRes(3).energyVal should approximate(-9.15) + + } withClue buildDebugString(batVars) + + model.release() + } + + "balance out additional power exceeding maximum battery power" in { + + given model: MPModel = MPModel(SolverLib.oJSolver) + + val constFlex = EnergyBoundariesFlexOptions( + AssetEnergyBoundaries( + Seq(5, -60, 110, -2).toPowerMap + ) + ) + + val flexOptions = Map( + loadUUID -> constFlex, + batUUID -> batFlex, + ) + + val assetVars = + OptimizedFlexStrat.addAssetConstraints(flexOptions, sampleTime, ticks) + + assetVars.toSeq should have length 2 + assetVars.foreach(_.results should have length 1) + assetVars.foreach(_.results.foreach(_ should have length 4)) + + val objectiveContainer = OptimizedFlexStrat.buildObjective( + assetVars, + zeroKW, + MinAbsPowerObjectiveFactory, + ) + + model.minimize(objectiveContainer.objective) + model.start(timeLimit = 10000) + + model.getStatus shouldBe SolutionStatus.OPTIMAL + + /* + EXPECTED RESULTS + Battery should be able to cover the additional power + up to its maximum power + */ + + val batVars = assetVars + .find(_.assetUuid == batUUID) + .getOrElse(fail(s"No asset variables for battery ($batUUID) found.")) + + val batRes = batVars.results.headOption + .getOrElse(fail(s"Empty results for battery ($batUUID).")) + + { + objectiveContainer.softConstraints.foreach { constraint => + withClue(constraint.getWarningMessage) { + constraint.getError should be < constraintTolerance + } + } + + // discharging 5 kWh plus 1.25 kWh losses + batRes(0).pVal should approximate(-5) + batRes(0).energyVal should approximate(-6.25) + + // charging 10 kWh minus 2 kWh losses + batRes(1).pVal should approximate(10) + batRes(1).energyVal should approximate(1.75) + + // discharging 10 kWh plus 2.5 kWh losses + batRes(2).pVal should approximate(-10) + batRes(2).energyVal should approximate(-10.75) + + // charging 2 kWh minus 0.4 kWh losses + batRes(3).pVal should approximate(2) + batRes(3).energyVal should approximate(-9.15) + + } withClue buildDebugString(batVars) + + model.release() + } + + "balance out additional power exceeding energy storage capacity" in { + + given model: MPModel = MPModel(SolverLib.oJSolver) + + val constFlex = EnergyBoundariesFlexOptions( + AssetEnergyBoundaries( + Seq(-10, -10, 10, 10).toPowerMap + ) + ) + + val flexOptions = Map( + loadUUID -> constFlex, + batUUID -> batFlex, + ) + + val assetVars = + OptimizedFlexStrat.addAssetConstraints(flexOptions, sampleTime, ticks) + + assetVars.toSeq should have length 2 + assetVars.foreach(_.results should have length 1) + assetVars.foreach(_.results.foreach(_ should have length 4)) + + val objectiveContainer = OptimizedFlexStrat.buildObjective( + assetVars, + zeroKW, + MinAbsPowerObjectiveFactory, + ) + + model.minimize(objectiveContainer.objective) + model.start(timeLimit = 10000) + + model.getStatus shouldBe SolutionStatus.OPTIMAL + + /* + EXPECTED RESULTS + Since excess power costs the same at all points in time and + at all magnitudes, there are many optimal solutions. + Thus, we only test for things that are true for every optimal + solution: We know when the battery should be definitely + full/empty and how much energy was charged/discharged. + + The soft constraints are vital here. Without them, + optimization would overestimate the losses in the first half + in order to achieve total power closer to zero. + */ + + val batVars = assetVars + .find(_.assetUuid == batUUID) + .getOrElse(fail(s"No asset variables for battery ($batUUID) found.")) + + val batRes = batVars.results.headOption + .getOrElse(fail(s"Empty results for battery ($batUUID).")) + + { + objectiveContainer.softConstraints.foreach { constraint => + withClue(constraint.getWarningMessage) { + constraint.getError should be < constraintTolerance + } + } + + // possibly charging + batRes(0).pVal should be >= 0d + batRes(0).energyVal should (be >= 0d and be <= 12d) + + // possibly charging, now we should have reached 24 kWh + // (12 kWh above starting energy) + batRes(1).pVal should be >= 0d + batRes(1).energyVal should approximate(12d) + + // we should've charged 12 kWh plus 3 kWh losses + batRes(0).pVal + batRes(1).pVal should approximate(15d) + + // possibly discharging + batRes(2).pVal should be <= 0d + batRes(2).energyVal should (be >= -12d and be <= 12d) + + // possibly discharging, now we should have reached 0 kWh + // (12 kWh below starting energy) + batRes(3).pVal should be <= 0d + batRes(3).energyVal should approximate(-12d) + + // we should've discharged 24 kWh minus 4.8 kWh losses + batRes(2).pVal + batRes(3).pVal should approximate(-19.2d) + + } withClue buildDebugString(batVars) + + model.release() + + } + + "balance out additional power exceeding maximum battery power and energy storage capacity" in { + + given model: MPModel = MPModel(SolverLib.oJSolver) + + val constFlex = EnergyBoundariesFlexOptions( + AssetEnergyBoundaries( + Seq(-10, -50, 20, 30).toPowerMap + ) + ) + + val flexOptions = Map( + loadUUID -> constFlex, + batUUID -> batFlex, + ) + + val assetVars = + OptimizedFlexStrat.addAssetConstraints(flexOptions, sampleTime, ticks) + + assetVars.toSeq should have length 2 + assetVars.foreach(_.results should have length 1) + assetVars.foreach(_.results.foreach(_ should have length 4)) + + val objectiveContainer = OptimizedFlexStrat.buildObjective( + assetVars, + zeroKW, + MinAbsPowerObjectiveFactory, + ) + + model.minimize(objectiveContainer.objective) + model.start(timeLimit = 10000) + + model.getStatus shouldBe SolutionStatus.OPTIMAL + + /* + EXPECTED RESULTS + Since excess power costs the same at all points in time and + at all magnitudes, there are many optimal solutions. + Thus, we only test for things that are true for every optimal + solution: We know when the battery should be definitely + full/empty and how much energy was charged/discharged. + + The soft constraints are vital here. Without them, + optimization would overestimate the losses in the first half + in order to achieve total power closer to zero. + */ + + val batVars = assetVars + .find(_.assetUuid == batUUID) + .getOrElse(fail(s"No asset variables for battery ($batUUID) found.")) + + val batRes = batVars.results.headOption + .getOrElse(fail(s"Empty results for battery ($batUUID).")) + + { + objectiveContainer.softConstraints.foreach { constraint => + withClue(constraint.getWarningMessage) { + constraint.getError should be < constraintTolerance + } + } + + // possibly charging + batRes(0).pVal should be >= 0d + batRes(0).energyVal should (be >= 0d and be <= 12d) + + // possibly charging, now we should have reached 24 kWh + // (12 kWh above starting energy) + batRes(1).pVal should be >= 0d + batRes(1).energyVal should approximate(12d) + + // we should've charged 12 kWh plus 3 kWh losses + batRes(0).pVal + batRes(1).pVal should approximate(15d) + + // possibly discharging + batRes(2).pVal should be <= 0d + batRes(2).energyVal should (be >= -12d and be <= 12d) + + // possibly discharging, now we should have reached 0 kWh + // (12 kWh below starting energy) + batRes(3).pVal should be <= 0d + batRes(3).energyVal should approximate(-12d) + + // we should've discharged 24 kWh minus 4.8 kWh losses + batRes(2).pVal + batRes(3).pVal should approximate(-19.2d) + + } withClue buildDebugString(batVars) + + model.release() + + } + + "balance out additional power exceeding energy storage capacity when discharging first" in { + + given model: MPModel = MPModel(SolverLib.oJSolver) + + val constFlex = EnergyBoundariesFlexOptions( + AssetEnergyBoundaries( + Seq(1, 1, -10, -10).toPowerMap + ) + ) + + val flexOptions = Map( + loadUUID -> constFlex, + batUUID -> batFlex, + ) + + val assetVars = + OptimizedFlexStrat.addAssetConstraints(flexOptions, sampleTime, ticks) + + assetVars.toSeq should have length 2 + assetVars.foreach(_.results should have length 1) + assetVars.foreach(_.results.foreach(_ should have length 4)) + + val objectiveContainer = OptimizedFlexStrat.buildObjective( + assetVars, + zeroKW, + MinAbsPowerObjectiveFactory, + ) + + model.minimize(objectiveContainer.objective) + model.start(timeLimit = 10000) + + model.getStatus shouldBe SolutionStatus.OPTIMAL + + /* + EXPECTED RESULTS + The soft constraints are vital here. Without them, + optimization would overestimate the losses of discharging + in the first half in order allow for more charging in the + second half. + */ + + val batVars = assetVars + .find(_.assetUuid == batUUID) + .getOrElse(fail(s"No asset variables for battery ($batUUID) found.")) + + val batRes = batVars.results.headOption + .getOrElse(fail(s"Empty results for battery ($batUUID).")) + + { + objectiveContainer.softConstraints.foreach { constraint => + withClue(constraint.getWarningMessage) { + constraint.getError should be < constraintTolerance + } + } + + // discharging 1 kWh plus 0.25 kWh losses + batRes(0).pVal should approximate(-1d) + batRes(0).energyVal should approximate(-1.25d) + + // discharging 1 kWh plus 0.25 kWh losses + batRes(1).pVal should approximate(-1d) + batRes(1).energyVal should approximate(-2.5d) + + // possibly charging + batRes(2).pVal should be >= 0d + batRes(2).energyVal should (be >= -2.5d and be <= 12d) + + // possibly charging, now we should have reached 24 kWh + // (12 kWh above starting energy) + batRes(3).pVal should be >= 0d + batRes(3).energyVal should approximate(12d) + + // we should've charged 14.5 kWh plus 3.625 kWh losses + batRes(2).pVal + batRes(3).pVal should approximate(18.125d) + + } withClue buildDebugString(batVars) + + model.release() + + } + + } + + "provided with battery flex options without losses" should { + + given ticks: Seq[Long] = + Range.Long.inclusive(0, 4 * sampleTicks, sampleTicks) + + // low efficiency for simplicity of the test + val batFlex = EnergyBoundariesFlexOptions( + AssetEnergyBoundaries( + eStorage = KilowattHours(24), + currentEnergy = KilowattHours(12), + pMax = Kilowatts(10), + etaCharge = Each(1), + etaDischarge = Each(1), + currentTick = 0L, + ) + ) + + // dummy energy conversion factor + given EnergyConversionFactor = + EnergyConversionFactor( + batFlex.energyBoundaries.headOption.value.etaCharge, + Each(1), + ) + + "balance out additional power within maximum battery power" in { + + given model: MPModel = MPModel(SolverLib.oJSolver) + + val constFlex = EnergyBoundariesFlexOptions( + AssetEnergyBoundaries( + Seq(5, -10, 10, -2).toPowerMap + ) + ) + + val flexOptions = Map( + loadUUID -> constFlex, + batUUID -> batFlex, + ) + + val assetVars = + OptimizedFlexStrat.addAssetConstraints(flexOptions, sampleTime, ticks) + + val objectiveContainer = OptimizedFlexStrat.buildObjective( + assetVars, + zeroKW, + MinAbsPowerObjectiveFactory, + ) + + model.minimize(objectiveContainer.objective) + model.start(timeLimit = 10000) + + model.getStatus shouldBe SolutionStatus.OPTIMAL + + /* + EXPECTED RESULTS + Battery should be able to fully cover the additional power. + No losses should be subtracted. + */ + + val batVars = assetVars + .find(_.assetUuid == batUUID) + .getOrElse(fail(s"No asset variables for battery ($batUUID) found.")) + + val batRes = batVars.results.headOption + .getOrElse(fail(s"Empty results for battery ($batUUID).")) + + { + objectiveContainer.softConstraints.foreach { constraint => + withClue(constraint.getWarningMessage) { + constraint.getError should be < constraintTolerance + } + } + + // discharging 5 kWh + batRes(0).pVal should approximate(-5) + batRes(0).energyVal should approximate(-5) + + // charging 10 kWh + batRes(1).pVal should approximate(10) + batRes(1).energyVal should approximate(5) + + // discharging 10 kWh + batRes(2).pVal should approximate(-10) + batRes(2).energyVal should approximate(-5) + + // charging 2 kWh + batRes(3).pVal should approximate(2) + batRes(3).energyVal should approximate(-3) + + } withClue buildDebugString(batVars) + + model.release() + } + + } + + "provided with energy boundary flex options and an objective factory" should { + + given ticks: Seq[Long] = + Range.Long.inclusive(0, 12 * sampleTicks, sampleTicks) + + // 33 kWh of feed-in in total, more than battery can store + val pvFlex = EnergyBoundariesFlexOptions( + AssetEnergyBoundaries( + Seq(0, -6, -8, -7, -12, 0, 0, 0, 0, 0, 0, 0).toPowerMap + ) + ) + + // 31 kWh of load in total, more than battery can provide + val loadFlex = EnergyBoundariesFlexOptions( + AssetEnergyBoundaries( + Seq(0, 0, 0, 0, 0, 0, 5, 10, 3, 7, 6, 0).toPowerMap + ) + ) + + // low efficiency for simplicity of the test + val batFlex = EnergyBoundariesFlexOptions( + AssetEnergyBoundaries( + eStorage = KilowattHours(20), + currentEnergy = KilowattHours(0), + pMax = Kilowatts(10), + etaCharge = Each(0.8), + etaDischarge = Each(0.8), + currentTick = 0L, + ) + ) + + val flexOptions = Map( + pvUUID -> pvFlex, + loadUUID -> loadFlex, + batUUID -> batFlex, + ) + + // since energy values have been adapted, we need this + // factor to convert back to "real" values + given EnergyConversionFactor = + EnergyConversionFactor( + batFlex.energyBoundaries.headOption.value.etaCharge, + OptimizedFlexStrat + .adaptEnergyBoundaries(batFlex.energyBoundaries.headOption.value) + .etaCharge, + ) + + "compensate fixed powers when using linear objective" in { + + given model: MPModel = MPModel(SolverLib.oJSolver) + + val assetVars = + OptimizedFlexStrat.addAssetConstraints(flexOptions, sampleTime, ticks) + + val objectiveContainer = OptimizedFlexStrat.buildObjective( + assetVars, + zeroKW, + MinAbsPowerObjectiveFactory, + ) + + model.minimize(objectiveContainer.objective) + model.start(timeLimit = 10000) + + model.getStatus shouldBe SolutionStatus.OPTIMAL + + /* + EXPECTED RESULTS + Since excess power costs the same at all points in time and + at all magnitudes, there are many optimal solutions. + Thus, we only test for things that are true for every optimal + solution: We know when the battery should be definitely + full/empty and how much energy was charged/discharged. + + The soft constraints are vital here. Without them, + optimization would overestimate the losses in the first half + in order to achieve total power closer to zero. + */ + + val batVars = assetVars + .find(_.assetUuid == batUUID) + .getOrElse(fail(s"No asset variables for battery ($batUUID) found.")) + + val batRes = batVars.results.headOption + .getOrElse(fail(s"Empty results for battery ($batUUID).")) + + { + objectiveContainer.softConstraints.foreach { constraint => + withClue(constraint.getWarningMessage) { + constraint.getError should be < constraintTolerance + } + } + + batRes should have size 12 + + batRes.slice(0, 4).foreach { + _.energyVal should (be >= 0d and be < 20d) + } + + batRes.slice(4, 6).foreach { + _.energyVal should approximate(20d) + } + + batRes.slice(6, 10).foreach { + _.energyVal should (be >= 0d and be < 20d) + } + + batRes.slice(10, 12).foreach { + _.energyVal should approximate(0d) + } + + // we should've charged 20 kWh plus 5 kWh losses + val inputCharged = batRes.slice(0, 6).map(_.pVal).sum + inputCharged should approximate(25) + + // we should've discharged 20 kWh minus 4 kWh losses + val outputDischarged = + batRes.slice(6, 12).map(_.pVal).sum + outputDischarged should approximate(-16d) + + } withClue buildDebugString(batVars) + + model.release() + + } + + "minimize peaks when using quadratic objective" in { + + given model: MPModel = MPModel(SolverLib.oJSolver) + + val assetVars = + OptimizedFlexStrat.addAssetConstraints(flexOptions, sampleTime, ticks) + + val objectiveContainer = OptimizedFlexStrat.buildObjective( + assetVars, + zeroKW, + LinearizedQuadraticPowerObjectiveFactory( + segmentCount = 10, + lastSegment = 10d, // 10 kW + ), + ) + + model.minimize(objectiveContainer.objective) + model.start(timeLimit = 10000) + + model.getStatus shouldBe SolutionStatus.OPTIMAL + + /* + EXPECTED RESULTS + When using the (linearized) quadratic objective, the + test is designed in a way such that only one optimal + solution exists: During the first phase (charging), + battery power is utilized so that 2 kW remains at every + time step. During the second phase (discharging), there's + exactly enough energy available to go down to 3 kW at + every time step. + + The soft constraints are vital here. Without them, + optimization would overestimate the losses in the first half + in order to achieve total power closer to zero. + */ + + val batVars = assetVars + .find(_.assetUuid == batUUID) + .getOrElse(fail(s"No asset variables for battery ($batUUID) found.")) + + val batRes = batVars.results.headOption + .getOrElse(fail(s"Empty results for battery ($batUUID).")) + + { + objectiveContainer.softConstraints.foreach { constraint => + withClue(constraint.getWarningMessage) { + constraint.getError should be < constraintTolerance + } + } + + // 0 kW to compensate + batRes(0).pVal should approximate(0d) + batRes(0).energyVal should approximate(0d) + + // 6 kW of feed-in to compensate, 2 kW remains + batRes(1).pVal should approximate(4d) + batRes(1).energyVal should approximate(3.2d) + + // 8 kW of feed-in to compensate, 2 kW remains + batRes(2).pVal should approximate(6d) + batRes(2).energyVal should approximate(8d) + + // 7 kW of feed-in to compensate, 2 kW remains + batRes(3).pVal should approximate(5d) + batRes(3).energyVal should approximate(12d) + + // 12 kW of feed-in to compensate, 2 kW remains + batRes(4).pVal should approximate(10d) + batRes(4).energyVal should approximate(20d) + + // 0 kW to compensate + batRes(5).pVal should approximate(0d) + batRes(5).energyVal should approximate(20d) + + // 5 kW of load to compensate, 3 kW remains + batRes(6).pVal should approximate(-2d) + batRes(6).energyVal should approximate(17.5d) + + // 10 kW of load to compensate, 3 kW remains + batRes(7).pVal should approximate(-7d) + batRes(7).energyVal should approximate(8.75d) + + // 3 kW of load to compensate, 3 kW remains + batRes(8).pVal should approximate(0d) + batRes(8).energyVal should approximate(8.75d) + + // 7 kW of load to compensate, 3 kW remains + batRes(9).pVal should approximate(-4d) + batRes(9).energyVal should approximate(3.75d) + + // 6 kW of load to compensate, 3 kW remains + batRes(10).pVal should approximate(-3d) + batRes(10).energyVal should approximate(0d) + + // 0 kW to compensate + batRes(11).pVal should approximate(0d) + batRes(11).energyVal should approximate(0d) + + // we should've charged 20 kWh plus 5 kWh losses + val inputCharged = batRes.slice(0, 6).map(_.pVal).sum + inputCharged should approximate(25) + + // we should've discharged 20 kWh minus 4 kWh losses + val outputDischarged = + batRes.slice(6, 12).map(_.pVal).sum + outputDischarged should approximate(-16d) + + } withClue buildDebugString(batVars) + + model.release() + + } + } + } + +} + +object OptimizedFlexStratSpec extends OptionValues { + + trait BatteryTesting { + val currentEnergy: Energy + + def charge(power: Power, duration: Time): BatteryTesting + } + + final case class ClassicModel( + override val currentEnergy: Energy, + eStorage: Energy, + pMax: Power, + etaCharging: Dimensionless, + etaDischarging: Dimensionless, + ) extends BatteryTesting { + def charge(power: Power, duration: Time): BatteryTesting = { + val netPower = + if power > zeroKW then power * etaCharging.toEach + else power * 1 / etaDischarging.toEach + + copy(currentEnergy = currentEnergy + netPower * duration) + } + } + + final case class AdaptedModel( + override val currentEnergy: Energy, + eStorage: Energy, + pMax: Power, + etaAvg: Dimensionless, + ) extends BatteryTesting { + def charge(power: Power, duration: Time): BatteryTesting = { + val newEnergy = + currentEnergy + (power - power.abs * (1 - etaAvg.toEach)) * duration + + copy(currentEnergy = newEnergy) + } + } + +} diff --git a/src/test/scala/edu/ie3/simona/model/participant/PvModelSpec.scala b/src/test/scala/edu/ie3/simona/model/participant/PvModelSpec.scala index 907b12d79f..038eb39da3 100644 --- a/src/test/scala/edu/ie3/simona/model/participant/PvModelSpec.scala +++ b/src/test/scala/edu/ie3/simona/model/participant/PvModelSpec.scala @@ -12,10 +12,14 @@ import edu.ie3.datamodel.models.input.system.characteristic.CosPhiFixed import edu.ie3.datamodel.models.input.{NodeInput, OperatorInput} import edu.ie3.datamodel.models.voltagelevels.GermanVoltageLevelUtils import edu.ie3.simona.model.participant.PvModel.{PvState, RadiationData} -import edu.ie3.simona.ontology.messages.flex.FlexType +import edu.ie3.simona.ontology.messages.flex.{ + EnergyBoundariesFlexOptions, + FlexType, +} import edu.ie3.simona.service.Data.SecondaryData.WeatherData import edu.ie3.simona.test.common.{DefaultTestData, UnitSpec, WeatherTestData} import edu.ie3.util.quantities.PowerSystemUnits.* +import edu.ie3.util.scala.quantities.DefaultQuantities.{zeroKW, zeroKWh} import edu.ie3.util.scala.quantities.{ ApparentPower, Kilovoltamperes, @@ -23,8 +27,9 @@ import edu.ie3.util.scala.quantities.{ } import org.locationtech.jts.geom.{Coordinate, GeometryFactory, Point} import org.scalatest.GivenWhenThen -import squants.Each -import squants.energy.{Power, Watts} +import squants.energy.{Power, WattHours, Watts} +import squants.time.Hours +import squants.{Each, Energy} import tech.units.indriya.quantity.Quantities.getQuantity import tech.units.indriya.unit.Units.* @@ -40,6 +45,7 @@ class PvModelSpec // testing tolerances private given Power = Watts(1e-6) + private given Energy = WattHours(1e-6) // build the NodeInputModel (which defines the location of the pv input model) // the NodeInputModel needs a GeoReference for the Pv to work @@ -177,15 +183,25 @@ class PvModelSpec val flexOptions = pvModel - .flexModels(FlexType.MathProgramming) + .flexModels(FlexType.EnergyBoundaries) .determineFlexOptions(state) flexOptions match { - case PowerSeriesMathFlexOptions(powerMap) => - powerMap should have size expectedResults.size - expectedResults.foreach { case (tick, expectedResult) => - powerMap(tick) should approximate(expectedResult) - } + case EnergyBoundariesFlexOptions(boundaries) => + boundaries should have size 1 + + val energyLimits = boundaries.headOption.value.energyLimits + energyLimits should have size expectedResults.size + 1 + + expectedResults + // adding dummy value so that last energy is tested + .appended(14400L -> zeroKW) + .foldLeft(zeroKWh) { case (expectedEnergy, (tick, expectedPower)) => + energyLimits(tick).getUpper should approximate(expectedEnergy) + energyLimits(tick).getLower should approximate(expectedEnergy) + + expectedEnergy + expectedPower * Hours(1) + } case unexpected => fail(s"Received unexpected flex options $unexpected") } diff --git a/src/test/scala/edu/ie3/simona/model/participant/WecModelSpec.scala b/src/test/scala/edu/ie3/simona/model/participant/WecModelSpec.scala index 50858c1637..fb1c70f85b 100644 --- a/src/test/scala/edu/ie3/simona/model/participant/WecModelSpec.scala +++ b/src/test/scala/edu/ie3/simona/model/participant/WecModelSpec.scala @@ -16,13 +16,18 @@ import edu.ie3.datamodel.models.input.system.characteristic.{ import edu.ie3.datamodel.models.input.{NodeInput, OperatorInput} import edu.ie3.datamodel.models.voltagelevels.GermanVoltageLevelUtils import edu.ie3.simona.model.participant.WecModel.{AirWeatherData, WecState} -import edu.ie3.simona.ontology.messages.flex.FlexType +import edu.ie3.simona.ontology.messages.flex.{ + EnergyBoundariesFlexOptions, + FlexType, +} import edu.ie3.simona.service.Data.SecondaryData.WeatherData import edu.ie3.simona.test.common.{UnitSpec, WeatherTestData} import edu.ie3.util.quantities.PowerSystemUnits -import squants.energy.{Power, Watts} +import edu.ie3.util.scala.quantities.DefaultQuantities.{zeroKW, zeroKWh} +import squants.energy.{Energy, Power, WattHours, Watts} import squants.motion.{MetersPerSecond, Pascals} import squants.thermal.Celsius +import squants.time.Hours import squants.{Dimensionless, Each} import tech.units.indriya.quantity.Quantities import tech.units.indriya.unit.Units.{METRE, PERCENT, SQUARE_METRE} @@ -32,9 +37,11 @@ import scala.collection.immutable.SortedMap class WecModelSpec extends UnitSpec with WeatherTestData { - protected given powerTolerance: Power = Watts(1e-6) - protected given dimensionlessTolerance: Dimensionless = Each(1e-12) - protected given doubleTolerance: Double = 1e-9 + // testing tolerances + private given Power = Watts(1e-6) + private given Energy = WattHours(1e-6) + private given Dimensionless = Each(1e-12) + private given Double = 1e-9 val nodeInput = new NodeInput( UUID.fromString("ad39d0b9-5ad6-4588-8d92-74c7d7de9ace"), @@ -204,15 +211,25 @@ class WecModelSpec extends UnitSpec with WeatherTestData { val flexOptions = wecModel - .flexModels(FlexType.MathProgramming) + .flexModels(FlexType.EnergyBoundaries) .determineFlexOptions(state) flexOptions match { - case PowerSeriesMathFlexOptions(powerMap) => - powerMap should have size expectedResults.size - expectedResults.foreach { case (tick, expectedResult) => - powerMap(tick) should approximate(expectedResult) - } + case EnergyBoundariesFlexOptions(boundaries) => + boundaries should have size 1 + + val energyLimits = boundaries.headOption.value.energyLimits + energyLimits should have size expectedResults.size + 1 + + expectedResults + // adding dummy value so that last energy is tested + .appended(43200L -> zeroKW) + .foldLeft(zeroKWh) { case (expectedEnergy, (tick, expectedPower)) => + energyLimits(tick).getUpper should approximate(expectedEnergy) + energyLimits(tick).getLower should approximate(expectedEnergy) + + expectedEnergy + expectedPower * Hours(1) + } case unexpected => fail(s"Received unexpected flex options $unexpected") } diff --git a/src/test/scala/edu/ie3/simona/model/participant/storage/StorageMathFlexModelSpec.scala b/src/test/scala/edu/ie3/simona/model/participant/storage/StorageMathFlexModelSpec.scala deleted file mode 100644 index 55b6fadd75..0000000000 --- a/src/test/scala/edu/ie3/simona/model/participant/storage/StorageMathFlexModelSpec.scala +++ /dev/null @@ -1,415 +0,0 @@ -/* - * © 2025. TU Dortmund University, - * Institute of Energy Systems, Energy Efficiency and Energy Economics, - * Research group Distribution grid planning and operation - */ - -package edu.ie3.simona.model.participant.storage - -import edu.ie3.simona.model.em.OptimizedFlexStrat -import edu.ie3.simona.model.participant.storage.StorageMathFlexModel.StorageMathFlexOptions -import edu.ie3.simona.model.participant.storage.StorageMathFlexModelSpec.* -import edu.ie3.simona.test.common.{MathFlexTestLike, UnitSpec} -import edu.ie3.util.scala.quantities.DefaultQuantities.{zeroKW, zeroKWh} -import optimus.algebra.{Double2Const, Zero} -import optimus.optimization.* -import optimus.optimization.enums.{SolutionStatus, SolverLib} -import optimus.optimization.model.MPFloatVar -import org.scalatest.OptionValues -import squants.energy.{Energy, KilowattHours, Kilowatts, WattHours} -import squants.time.Hours -import squants.{Dimensionless, Each, Power, Time} - -import java.util.UUID - -class StorageMathFlexModelSpec extends UnitSpec with MathFlexTestLike { - - // Testing tolerances - given Double = 1e-6 - given Energy = WattHours(1e-9) - val constraintTolerance = 1e-3 - - "StorageMathFlexModelSpec" should { - - "be created with correct parameters" in { - - val currentEnergy = KilowattHours(10) - val eStorage = KilowattHours(20) - val pMax = Kilowatts(10) - - val etas = Seq(.6, .65, .7, .75, .8, .85, 0.9, .92, .95, .98, 1) - - forEvery(Table("etaCharging", etas*)) { etaCharging => - forEvery(Table("etaDischarging", etas*)) { etaDischarging => - - val classic = ClassicModel( - currentEnergy = currentEnergy, - eStorage = eStorage, - pMax = pMax, - etaCharging = Each(etaCharging), - etaDischarging = Each(etaDischarging), - ) - - // charging until full, with maximum power - val power1 = pMax - val duration1 = (eStorage - currentEnergy) / (power1 * etaCharging) - // discharging three quarters, with half power - val power2 = -pMax / 2 - val time2 = eStorage * 0.75 / (-power2 * 1 / etaDischarging) - // charging until half, with quarter power - val power3 = pMax / 4 - val time3 = eStorage / 4 / (power3 * etaCharging) - // discharging until empty, with maximum power - val power4 = -pMax - val time4 = eStorage / 2 / (-power4 * 1 / etaDischarging) - - val adaptedFlexOptions = - StorageMathFlexOptions.createAdaptedFlexOptions( - currentEnergy, - eStorage, - pMax, - Each(etaCharging), - Each(etaDischarging), - ) - - val adapted = AdaptedModel( - currentEnergy = adaptedFlexOptions.currentEnergy, - eStorage = adaptedFlexOptions.eStorage, - pMax = adaptedFlexOptions.pMax, - etaAvg = adaptedFlexOptions.eta, - ) - - val results = IndexedSeq(classic, adapted).map( - _.charge(power1, duration1) - .charge(power2, time2) - .charge(power3, time3) - .charge(power4, time4) - .currentEnergy - ) - - // battery should be exactly empty - results(0) should approximate(zeroKWh) - results(1) should approximate(zeroKWh) - - } - } - } - - "balance out additional power with zero excess" in { - - // low efficiency for simplicity of the test - val fo = StorageMathFlexOptions.createAdaptedFlexOptions( - currentEnergy = KilowattHours(50), - eStorage = KilowattHours(100), - pMax = Kilowatts(10), - etaCharging = Each(0.8), - etaDischarging = Each(0.8), - ) - - val stepResolution = Hours(1) - val tickResolution = stepResolution.toSeconds.toLong - - given model: MPModel = MPModel(SolverLib.oJSolver) - - // since energy values have been adapted, we need this - // factor to convert back to "real" values - given EnergyConversionFactor = EnergyConversionFactor(Each(0.8), fo.eta) - - val container = OptimizedFlexStrat.addAssetConstraints( - assetUuid = UUID.randomUUID(), - flexOptions = fo, - ticks = Range.Long(0, tickResolution * 5, tickResolution), - ) - - container.states should have length 5 - container.operationVars should have length 4 - - // additional powers for each time step, all within pMax - val addPower = Seq(5d, -10d, 10d, -2d) - - val mainObjectiveDifferences = - container.operationVars.zip(addPower).map { case (opVar, add) => - val d = MPFloatVar.positive("d") - model.add(d >:= opVar.getPowerExpression + add) - model.add(d >:= -(opVar.getPowerExpression + add)) - d - } - - val softConstraints = - container.operationVars.flatMap(_.getSoftConstraint(stepResolution)) - - val objective = mainObjectiveDifferences - .appendedAll(softConstraints.map(_.getExpression)) - .reduceLeftOption(_ + _) - .getOrElse(Zero) - - model.minimize(objective) - - model.start() - - model.getStatus shouldBe SolutionStatus.OPTIMAL - - // Battery should be able to fully cover the additional power - - { - softConstraints.foreach { constraint => - withClue(constraint.getWarningMessage) { - constraint.getError should be < constraintTolerance - } - } - - container.states(0).energyVal should approximate(50) - - // discharging 5 kWh plus 1.25 kWh losses - container.operationVars(0).pVal should approximate(-5) - container.states(1).energyVal should approximate(43.75) - - // charging 10 kWh minus 2 kWh losses - container.operationVars(1).pVal should approximate(10) - container.states(2).energyVal should approximate(51.75) - - // discharging 10 kWh plus 2.5 kWh losses - container.operationVars(2).pVal should approximate(-10) - container.states(3).energyVal should approximate(39.25) - - // charging 2 kWh minus 0.4 kWh losses - container.operationVars(3).pVal should approximate(2) - container.states(4).energyVal should approximate(40.85) - - } withClue buildDebugString(container) - - model.release() - } - - "balance out additional power with large excess" in { - - // low efficiency for simplicity of the test - val fo = StorageMathFlexOptions.createAdaptedFlexOptions( - currentEnergy = KilowattHours(50), - eStorage = KilowattHours(100), - pMax = Kilowatts(10), - etaCharging = Each(0.8), - etaDischarging = Each(0.8), - ) - - val stepResolution = Hours(1) - val tickResolution = stepResolution.toSeconds.toLong - - given model: MPModel = MPModel(SolverLib.oJSolver) - - // since energy values have been adapted, we need this - // factor to convert back to "real" values - given EnergyConversionFactor = EnergyConversionFactor(Each(0.8), fo.eta) - - val container = OptimizedFlexStrat.addAssetConstraints( - assetUuid = UUID.randomUUID(), - flexOptions = fo, - ticks = Range.Long(0, tickResolution * 5, tickResolution), - ) - - container.states should have length 5 - container.operationVars should have length 4 - - // additional powers for each time step, some far beyond pMax - val addPower = Seq(5d, -60d, 110d, -2d) - - val mainObjectiveDifferences = - container.operationVars.zip(addPower).map { case (opVar, add) => - val d = MPFloatVar.positive("d") - model.add(d >:= opVar.getPowerExpression + add) - model.add(d >:= -(opVar.getPowerExpression + add)) - d - } - - val softConstraints = - container.operationVars.flatMap(_.getSoftConstraint(stepResolution)) - - val objective = mainObjectiveDifferences - .appendedAll(softConstraints.map(_.getExpression)) - .reduceLeftOption(_ + _) - .getOrElse(Zero) - - model.minimize(objective) - - model.start() - - model.getStatus shouldBe SolutionStatus.OPTIMAL - - // Battery should be able to fully cover the additional power - - { - softConstraints.foreach { constraint => - withClue(constraint.getWarningMessage) { - constraint.getError should be < constraintTolerance - } - } - - container.states(0).energyVal should approximate(50) - - // discharging 5 kWh plus 1.25 kWh losses - container.operationVars(0).pVal should approximate(-5) - container.states(1).energyVal should approximate(43.75) - - // charging 10 kWh minus 2 kWh losses - container.operationVars(1).pVal should approximate(10) - container.states(2).energyVal should approximate(51.75) - - // discharging 10 kWh plus 2.5 kWh losses - container.operationVars(2).pVal should approximate(-10) - container.states(3).energyVal should approximate(39.25) - - // charging 2 kWh minus 0.4 kWh losses - container.operationVars(3).pVal should approximate(2) - container.states(4).energyVal should approximate(40.85) - - } withClue buildDebugString(container) - - model.release() - } - - "work correctly at extreme values" in { - // low efficiency for simplicity of the test - val fo = StorageMathFlexOptions.createAdaptedFlexOptions( - currentEnergy = KilowattHours(10), - eStorage = KilowattHours(20), - pMax = Kilowatts(10), - etaCharging = Each(0.8), - etaDischarging = Each(0.8), - ) - - val stepResolution = Hours(1) - val tickResolution = stepResolution.toSeconds.toLong - - given model: MPModel = MPModel(SolverLib.oJSolver) - - // since energy values have been adapted, we need this - // factor to convert back to "real" values - given EnergyConversionFactor = EnergyConversionFactor(Each(0.8), fo.eta) - - val container = OptimizedFlexStrat.addAssetConstraints( - assetUuid = UUID.randomUUID(), - flexOptions = fo, - ticks = Range.Long(0, tickResolution * 5, tickResolution), - ) - - container.states should have length 5 - container.operationVars should have length 4 - - // additional powers for each time step - val addPower = Seq(-5d, -10d, 10d, 10d) - - val mainObjectiveDifferences = - container.operationVars.zip(addPower).map { case (opVar, add) => - val d = MPFloatVar.positive("d") - model.add(d >:= opVar.getPowerExpression + add) - model.add(d >:= -(opVar.getPowerExpression + add)) - d - } - - val softConstraints = - container.operationVars.flatMap(_.getSoftConstraint(stepResolution)) - - val objective = mainObjectiveDifferences - .appendedAll(softConstraints.map(_.getExpression)) - .reduceLeftOption(_ + _) - .getOrElse(Zero) - - model.minimize(objective) - - model.start() - - model.getStatus shouldBe SolutionStatus.OPTIMAL - - /* - EXPECTED RESULTS - Since excess power costs the same at all points in time and - at all magnitudes, there are many optimal solutions. - Thus, we only test for things that are true for every optimal - solution: We know when the battery should be definitely - full/empty and how much energy was charged/discharged. - */ - - { - softConstraints.foreach { constraint => - withClue(constraint.getWarningMessage) { - constraint.getError should be < constraintTolerance - } - } - - container.states(0).energyVal should approximate(10) - - // possibly charging - container.operationVars(0).pVal should be >= 0d - container.states(1).energyVal should (be >= 0d and be <= 20d) - - // possibly charging, now we should have reached 20 kWh - container.operationVars(1).pVal should be >= 0d - container.states(2).energyVal should approximate(20d) - - // we should've charged 10 kWh plus 2.5 kWh losses - val totalCharged = - container.operationVars(0).pVal + container.operationVars(1).pVal - totalCharged should approximate(12.5d) - - // possibly discharging - container.operationVars(2).pVal should be <= 0d - container.states(3).energyVal should (be >= 0d and be <= 20d) - - // possibly discharging, now we should have reached 0 kWh - container.operationVars(3).pVal should be <= 0d - container.states(4).energyVal should approximate(0d) - - // we should've discharged 20 kWh minus 4 kWh losses - val totalDischarged = - container.operationVars(2).pVal + container.operationVars(3).pVal - totalDischarged should approximate(-16d) - - } withClue buildDebugString(container) - - model.release() - - } - - } - -} - -object StorageMathFlexModelSpec extends OptionValues { - - trait BatteryTesting { - val currentEnergy: Energy - - def charge(power: Power, duration: Time): BatteryTesting - } - - final case class ClassicModel( - override val currentEnergy: Energy, - eStorage: Energy, - pMax: Power, - etaCharging: Dimensionless, - etaDischarging: Dimensionless, - ) extends BatteryTesting { - def charge(power: Power, duration: Time): BatteryTesting = { - val netPower = - if power > zeroKW then power * etaCharging.toEach - else power * 1 / etaDischarging.toEach - - copy(currentEnergy = currentEnergy + netPower * duration) - } - } - - final case class AdaptedModel( - override val currentEnergy: Energy, - eStorage: Energy, - pMax: Power, - etaAvg: Dimensionless, - ) extends BatteryTesting { - def charge(power: Power, duration: Time): BatteryTesting = { - val newEnergy = - currentEnergy + (power - power.abs * (1 - etaAvg.toEach)) * duration - - copy(currentEnergy = newEnergy) - } - } - -} diff --git a/src/test/scala/edu/ie3/simona/test/common/MathFlexTestLike.scala b/src/test/scala/edu/ie3/simona/test/common/OptimizingTestLike.scala similarity index 53% rename from src/test/scala/edu/ie3/simona/test/common/MathFlexTestLike.scala rename to src/test/scala/edu/ie3/simona/test/common/OptimizingTestLike.scala index c2842ed482..7b8557c9ae 100644 --- a/src/test/scala/edu/ie3/simona/test/common/MathFlexTestLike.scala +++ b/src/test/scala/edu/ie3/simona/test/common/OptimizingTestLike.scala @@ -6,37 +6,40 @@ package edu.ie3.simona.test.common -import edu.ie3.simona.model.em.OptimizedFlexStrat.AssetVarContainer -import edu.ie3.simona.model.participant.storage.StorageMathFlexModel.{ - StorageOperationVars, - StorageStateVars, +import edu.ie3.simona.model.em.opt.OptimizedFlexStrat.{ + AssetVarContainer, + StepResults, } import optimus.algebra.Const import optimus.optimization.model.MPVar +import org.scalatest.Assertions.fail import org.scalatest.OptionValues.convertOptionToValuable import squants.{Dimensionless, Power} -import squants.energy.Kilowatts +import squants.energy.{KilowattHours, Kilowatts} import scala.collection.immutable.SortedMap -trait MathFlexTestLike { +trait OptimizingTestLike { extension (seq: Seq[Int]) def toPowerMap(using ticks: Seq[Long]): SortedMap[Long, Power] = SortedMap.from(ticks.zip(seq.map(Kilowatts.apply))) - extension (state: StorageStateVars) + extension (res: StepResults) def energyVal(using conversion: EnergyConversionFactor): Double = { - val solution = state.storedEnergy match { - case constant: Const => constant.value - case mpVar: MPVar => mpVar.value.value + val solution = res.state + .getOrElse(fail("No state provided in StepResults!")) match { + case variable: MPVar => + variable.value.value + case const: Const => + const.value } + solution * conversion.factor } - extension (state: StorageOperationVars) def pVal: Double = - state.p.value.value + res.getOperationResult.toKilowatts final case class EnergyConversionFactor(factor: Double) @@ -49,10 +52,18 @@ trait MathFlexTestLike { } def buildDebugString( - batVars: AssetVarContainer[StorageStateVars, StorageOperationVars] + assetVars: AssetVarContainer )(using EnergyConversionFactor): String = - "\n\tDEBUGGING:" + - s"\n\t\tOperation values: ${batVars.operationVars.map(_.pVal).mkString(", ")}" + - s"\n\t\tState values: ${batVars.states.map(_.energyVal).mkString(", ")}" + s"\n\tDEBUGGING asset ${assetVars.assetUuid}:" + + assetVars.results + .map { res => + s"\n\t\tTrajectory: ${res + .map(step => + step.getOperationResult.toString + + step.state.map(_ => s" ( -> ${KilowattHours(step.energyVal).toString})").getOrElse("") + ) + .mkString(", ")}" + } + .mkString("") }