Skip to content

Commit

Permalink
pumping: Add pre-conditions for states, enforce more coding conventions
Browse files Browse the repository at this point in the history
  • Loading branch information
wasowski committed Nov 29, 2023
1 parent 9dcef9e commit 7a06e37
Show file tree
Hide file tree
Showing 3 changed files with 75 additions and 63 deletions.
110 changes: 61 additions & 49 deletions src/main/scala/symsim/examples/concrete/pumping/Pump.scala
Original file line number Diff line number Diff line change
@@ -1,22 +1,17 @@
package symsim
package examples.concrete.pumping
/** A case study of control of a water extraction at Marbjerg (operated by HOFOR).
*
* Source: Andreas Holck Høeg-Petersen. Reinforcement Learning for Controlling
* Groundwater Extraction. MSc thesis at the IT University of Copenhagen.
* August 2021. (esp. Chapters 2 and 3)
*/

import symsim.concrete.Randomized
import scala.math


/** Constants as fit in the model of AHHP (reference above).
*
* They are called a, b, c, and d parameters in his report.
*/
val c1: Double = 1.0 / 12
val c2: Double = 0.15 / (c1*80*12*24*365)
val c3: Double = 0.1 / (12*24*365)
val c2: Double = 0.15 / (c1 * 80 * 12 * 24 * 365)
val c3: Double = 0.1 / (12 * 24 * 365)
val c4: Double = 0.05

/** Size of the history window for the head-level */
Expand All @@ -33,6 +28,14 @@ case class PumpState (
w: Double, // Water level in the ground deposits (unobservable)
phm: List[Double]): // Past head means (a sliding window of history)

require (f >= FlowMin, s"Flow is too low: $f")
require (f <= FlowMax, s"Flow is too high: $f")
require (h <= HeadMax, s"Head is too high: $h")
require (hm >= HeadMin, s"Head mean is too low: $hm")
require (hm <= HeadMax, s"Head mean is too high: $hm")
require (w >= WaterMin, s"Water level is too low: $w")
require (w <= WaterMax, s"Water level is too high: $w")

override def toString: String =
s"[flow=$f, head=$h, head mean=$hm, tank level=$tl, time=$t, "
+ s"water=$w, past head means=$phm]"
Expand All @@ -42,9 +45,9 @@ end PumpState

/** Discretize cut points for flow, head, and tank variables./
*
* We observe to the first point on the list that is lower than the actual
* value of the variable. So we 'round down'. Values below the last element
* in the cut-point list are rounded up (to the last value).
* We discretize to the first point on the list that is lower than the actual
* value of the variable. So we 'round down'. Values below the last element in
* the cut-point list are rounded up (to the last value).
*
* For now, we have preconditions about low values in discretization so the
* getOrElse (rounding up) should never happen.
Expand All @@ -57,24 +60,24 @@ def closest (value: Double) (cutPoints: List[Double]): Double =
.find { value >= _ }
.getOrElse (cutPoints.last)

val HEAD_MIN: Double = 7.0
val HEAD_MAX: Double = 10.84
val FLOW_MIN: Double = 0.0
val FLOW_MAX: Double = 120.0
val TANK_MIN: Double = 0.0
val TANK_MAX: Double = 2000.0
val WATER_MIN: Double = 0.0
val WATER_MAX: Double = 9.11
val HeadMin = 7.0
val HeadMax = 10.84
val FlowMin = 0.0
val FlowMax = 120.0
val TankMin = 0.0
val TankMax = 2000.0
val WaterMin = 0.0
val WaterMax = 9.11

val flowCutPoints = List (FLOW_MAX, 115.0, 110.0, 105.0, 100.0, 95.0,
90.0, 85.0, 80.0, 75.0, 70.0, 65.0, 60.0, 55.0, 50.0, FLOW_MIN)
val FlowCutPoints = List (FlowMax, 115.0, 110.0, 105.0, 100.0, 95.0,
90.0, 85.0, 80.0, 75.0, 70.0, 65.0, 60.0, 55.0, 50.0, FlowMin)

val headCutPoints = List (HEAD_MAX, 10.68, 10.52, 10.36, 10.2, 10.04,
val HeadCutPoints = List (HeadMax, 10.68, 10.52, 10.36, 10.2, 10.04,
9.88, 9.72, 9.56, 9.4, 9.24, 9.08, 8.92, 8.76, 8.6, 8.44, 8.28,
8.12, 7.96, 7.8, 7.64, 7.48, 7.32, 7.16, HEAD_MIN)
8.12, 7.96, 7.8, 7.64, 7.48, 7.32, 7.16, HeadMin)

val tankCutPoints = List (TANK_MAX, 1800.0, 1600.0, 1400.0,
1200.0, 1000.0, 800.0, 600.0, 400.0, 200.0, TANK_MIN)
val TankCutPoints = List (TankMax, 1800.0, 1600.0, 1400.0,
1200.0, 1000.0, 800.0, 600.0, 400.0, 200.0, TankMin)

/** The (discrete) state, observable by the pump controller. */
case class ObservablePumpState (
Expand All @@ -83,6 +86,11 @@ case class ObservablePumpState (
hm: Double, // Discretized, derived head mean
tl: Double): // Discretized observable tank level

require (FlowCutPoints.contains (f), s"Wrong discrete flow: $f")
require (HeadCutPoints.contains (h), s"Wrong discrete head level: $h")
require (HeadCutPoints.contains (hm), s"Wrong discrete mean head level: $hm")
require (TankCutPoints.contains (tl), s"Wrong discrete tank level: $tl")

override def toString: String =
s"(flow=$f, head=$h, head mean=$hm, tank level=$tl)"

Expand All @@ -93,6 +101,12 @@ type PumpAction = Double
type PumpReward = Double


/** A case study of control of a water extraction at Marbjerg (operated by HOFOR).
*
* Source: Andreas Holck Høeg-Petersen. Reinforcement Learning for Controlling
* Groundwater Extraction. MSc thesis at the IT University of Copenhagen.
* August 2021. (esp. Chapters 2 and 3)
*/
object Pump extends
Agent[PumpState, ObservablePumpState, PumpAction, PumpReward, Randomized],
Episodic:
Expand All @@ -103,32 +117,30 @@ object Pump extends

/** If isFinal is true in a state, then the episode is complete */
def isFinal (s: PumpState): Boolean =
s.t >= 4000 || s.h < HEAD_MIN || s.tl > TANK_MAX || s.tl < TANK_MIN
s.t >= 4000 || s.h < HeadMin || s.tl > TankMax || s.tl < TankMin


def observe (s: PumpState): ObservablePumpState =
require (s.tl >= TANK_MIN, s"s.tl = ${s.tl} >= $TANK_MIN")

val df = closest (s.f) (flowCutPoints)
val dh = closest (s.h) (headCutPoints)
val dhm = closest (s.hm) (headCutPoints)
val dtl = closest (s.tl) (tankCutPoints)

require (s.tl >= TankMin, s"s.tl = ${s.tl} >= $TankMin")
val df = closest (s.f) (FlowCutPoints)
val dh = closest (s.h) (HeadCutPoints)
val dhm = closest (s.hm) (HeadCutPoints)
val dtl = closest (s.tl) (TankCutPoints)
ObservablePumpState (df, dh, dhm, dtl)


def reward (source: PumpState) (target: PumpState) (a: PumpAction): Double =
if target.h < HEAD_MIN || target.tl < TANK_MIN || target.tl > TANK_MAX
if target.h < HeadMin || target.tl < TankMin || target.tl > TankMax
then -9999
else
val flowReward = if target.f != source.f then -0.5 else 0
val headCost = (1 + (target.h - target.hm).abs) * (1 + (target.h - target.hm).abs)
flowReward - headCost


val amp: Double = 0.41852857594808646
val freq: Double = 0.0000597030105413
val phase: Double = -6347.109214682171
val Amp = 0.41852857594808646
val Freq = 0.0000597030105413
val Phase = -6347.109214682171


override def step (s: PumpState) (a: PumpAction)
Expand All @@ -141,10 +153,10 @@ object Pump extends
cd <- getDemand (s.t%24 + 1)
d = cd + nd
tl1 = s.tl + c1 * (f1 - d)
h1 = s.h + c4 * (s.w + (c1*f1 / Math.PI))
h1 = s.h + c4 * (s.w + (c1*f1 / math.Pi))
nw <- Randomized.gaussian (0.0, 1.0)
w1 = s.w - c2*(c1*f1) + c3 +
(amp*Math.sin (2*Math.PI*(s.t + phase) / freq)) + nw
(Amp * math.sin (2 * math.Pi*(s.t + Phase) / Freq)) + nw
hm1 = (1.0 / k)*s.phm.sum
phm1 = (s.hm:: s.phm).slice (0, k)
s1 = PumpState (f1, h1, hm1, tl1, s.t + 1, w1, phm1)
Expand Down Expand Up @@ -196,10 +208,10 @@ object PumpInstances

given enumState: BoundedEnumerable[ObservablePumpState] =
val ss: List[ObservablePumpState] = for
f <- flowCutPoints
h <- headCutPoints
hm <- headCutPoints
tl <- tankCutPoints
f <- FlowCutPoints
h <- HeadCutPoints
hm <- HeadCutPoints
tl <- TankCutPoints
yield ObservablePumpState (f, h, hm, tl)
BoundedEnumerableFromList (ss*)

Expand All @@ -213,13 +225,13 @@ object PumpInstances
concrete.Randomized.canTestInRandomized

lazy val genPumpState: Gen[PumpState] = for
f <- Gen.choose[Double] (FLOW_MIN, FLOW_MAX)
h <- Gen.choose[Double] (HEAD_MIN, HEAD_MAX)
hm <- Gen.choose[Double] (HEAD_MIN, HEAD_MAX)
tl <- Gen.choose[Double] (TANK_MIN, TANK_MAX)
f <- Gen.choose[Double] (FlowMin, FlowMax)
h <- Gen.choose[Double] (HeadMin, HeadMax)
hm <- Gen.choose[Double] (HeadMin, HeadMax)
tl <- Gen.choose[Double] (TankMin, TankMax)
t <- Gen.choose[Int] (0, 24)
w <- Gen.choose[Double] (WATER_MIN, WATER_MAX)
phm <- Gen.listOfN (5, Gen.choose (HEAD_MIN, HEAD_MAX))
w <- Gen.choose[Double] (WaterMin, WaterMax)
phm <- Gen.listOfN (5, Gen.choose (HeadMin, HeadMax))
yield PumpState (f, h, hm, tl, t, w, phm)

given arbitraryState: Arbitrary[PumpState] = Arbitrary (genPumpState)
Expand Down
Original file line number Diff line number Diff line change
@@ -1,8 +1,8 @@
package symsim
package examples.concrete.pumping

import laws.AgentLaws
import laws.EpisodicLaws
import symsim.laws.AgentLaws
import symsim.laws.EpisodicLaws

class PumpIsAgentSpec
extends SymSimSpec:
Expand Down
24 changes: 12 additions & 12 deletions src/test/scala/symsim/examples/concrete/pumping/PumpSpec.scala
Original file line number Diff line number Diff line change
Expand Up @@ -2,8 +2,8 @@ package symsim
package examples.concrete.pumping

import symsim.concrete.Randomized.given
import CanTestIn.given
import CanTestIn.*
import symsim.CanTestIn.given
import symsim.CanTestIn.*

import org.scalacheck.{Arbitrary, Gen, Prop}
import org.scalacheck.Arbitrary.*
Expand All @@ -19,13 +19,13 @@ class PumpSpec
extends org.scalacheck.Properties ("Pumping"):

// Generators of test data
val flow = Gen.choose[Double] (FLOW_MIN, FLOW_MAX)
val head = Gen.choose[Double] (HEAD_MIN, HEAD_MAX)
val head_mean = Gen.choose[Double] (HEAD_MIN, HEAD_MAX)
val tank = Gen.choose[Double] (TANK_MIN, TANK_MAX)
val flow = Gen.choose[Double] (FlowMin, FlowMax)
val head = Gen.choose[Double] (HeadMin, HeadMax)
val head_mean = Gen.choose[Double] (HeadMin, HeadMax)
val tank = Gen.choose[Double] (TankMin, TankMax)
val time = Gen.choose[Int] (0, 24)
val water = Gen.choose[Double] (WATER_MIN, WATER_MAX)
val past_head_mean = Gen.listOfN (5, Gen.choose (HEAD_MIN, HEAD_MAX))
val water = Gen.choose[Double] (WaterMin, WaterMax)
val past_head_mean = Gen.listOfN (5, Gen.choose (HeadMin, HeadMax))
val actions = Gen.oneOf (Pump.instances.enumAction.membersAscending)
val obsStates = Gen.oneOf (Pump.instances.allObservableStates)

Expand All @@ -36,7 +36,7 @@ class PumpSpec
forAll (flow, head, head_mean, tank, time, water, past_head_mean, actions) {
(f, h, hm, tl, t, w, phm, a) =>
for (s1, r) <- Pump.step (PumpState (f, h, hm, tl, t, w, phm)) (a)
yield s1.tl >= TANK_MIN
yield s1.tl >= TankMin
}
}

Expand All @@ -48,9 +48,9 @@ class PumpSpec
}
}

property ("The time can be more than TANK_MAX (comparing is not correct)") = {
property ("The time can be more than TankMax (comparing is not correct)") = {
exists (time) { t =>
t > TANK_MAX
t > TankMax
}
}

Expand All @@ -65,7 +65,7 @@ class PumpSpec
(f, h, hm, tl, t, w, phm) =>
exists(actions) { a =>
for (s1, r) <- Pump.step(PumpState(f, h, hm, tl, t, w, phm))(a)
yield (s1.tl <= TANK_MAX)
yield (s1.tl <= TankMax)
}
}
}

0 comments on commit 7a06e37

Please sign in to comment.