diff --git a/src/main/scala/symsim/examples/concrete/pumping/Pump.scala b/src/main/scala/symsim/examples/concrete/pumping/Pump.scala index c9a2b072..e76c5c4a 100644 --- a/src/main/scala/symsim/examples/concrete/pumping/Pump.scala +++ b/src/main/scala/symsim/examples/concrete/pumping/Pump.scala @@ -1,13 +1,8 @@ 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). @@ -15,12 +10,12 @@ import symsim.concrete.Randomized * 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 */ -val k: Int = 5 +val MemoryWindowSize: Int = 5 /** Complete state of a pump, both observable and unobservable. */ @@ -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]" @@ -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. @@ -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 ( @@ -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)" @@ -93,70 +101,66 @@ 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: - - /** An upper bound on episode duration. Used only in testing */ - val TimeHorizon: Int = 20000 - + Agent[PumpState, ObservablePumpState, PumpAction, PumpReward, Randomized]: /** 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 - then -9999 + if isFinal (target) then -9999.0 else - val flowReward = if target.f != source.f then -0.5 else 0 + val flowReward = if target.f != source.f then -0.5 else 0.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) : Randomized[(PumpState, PumpReward)] = require (instances.enumAction.membersAscending.contains (a)) for - nf <- Randomized.gaussian (0.0, 1.0) - f1 = a + nf - nd <- Randomized.gaussian (0.1, 0.01) - cd <- getDemand (s.t%24 + 1) - d = cd + nd + f1 <- Randomized.gaussian (a, 1.0) + cd <- demand (s.t % 24 + 1) + d <- Randomized.gaussian (cd + 0.1, 0.01) tl1 = s.tl + c1 * (f1 - d) - 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 - hm1 = (1.0 / k)*s.phm.sum - phm1 = (s.hm:: s.phm).slice (0, k) + 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 + hm1 = (1.0/MemoryWindowSize) * s.phm.sum + phm1 = (s.hm:: s.phm).slice (0, MemoryWindowSize) s1 = PumpState (f1, h1, hm1, tl1, s.t + 1, w1, phm1) pr = reward (s) (s1) (a) yield (s1, pr) - def getDemand (t: Int): Randomized[Double] = + def demand (t: Int): Randomized[Double] = require (t >= 0 && t <= 24) if t < 5 then Randomized.between (5.0, 15.0) - if t < 12 then Randomized.between (15.0, 45.0) - if t < 22 then Randomized.between (20.0, 38.0) + else if t < 12 then Randomized.between (15.0, 45.0) + else if t < 22 then Randomized.between (20.0, 38.0) else Randomized.between (5.0, 20.0) @@ -196,10 +200,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*) @@ -213,13 +217,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) diff --git a/src/main/scala/symsim/laws/AgentLaws.scala b/src/main/scala/symsim/laws/AgentLaws.scala index a47aaa5a..2b7b720e 100644 --- a/src/main/scala/symsim/laws/AgentLaws.scala +++ b/src/main/scala/symsim/laws/AgentLaws.scala @@ -22,7 +22,8 @@ case class AgentLaws[State, ObservableState, Action, Reward, Scheduler[_]] import agent.instances.given - val observableStates = agent.instances.enumState.membersAscending + val observableStates = + agent.instances.enumState.membersAscending.toSet val laws: RuleSet = SimpleRuleSet ( "agent", diff --git a/src/test/scala/symsim/examples/concrete/pumping/Experiments.scala b/src/test/scala/symsim/examples/concrete/pumping/Experiments.scala index 74af5ecc..e79a7d38 100644 --- a/src/test/scala/symsim/examples/concrete/pumping/Experiments.scala +++ b/src/test/scala/symsim/examples/concrete/pumping/Experiments.scala @@ -2,10 +2,9 @@ package symsim package examples.concrete.pumping class Experiments - extends ExperimentSpec[PumpState,ObservablePumpState,PumpAction]: + extends ExperimentSpec[PumpState, ObservablePumpState, PumpAction]: // Import evidence that states and actions can be enumerated - import Pump.* import Pump.instances.given val sarsa = symsim.concrete.ConcreteSarsa ( diff --git a/src/test/scala/symsim/examples/concrete/pumping/PumpIsAgentSpec.scala b/src/test/scala/symsim/examples/concrete/pumping/PumpIsAgentSpec.scala index ed68f84c..32356c33 100644 --- a/src/test/scala/symsim/examples/concrete/pumping/PumpIsAgentSpec.scala +++ b/src/test/scala/symsim/examples/concrete/pumping/PumpIsAgentSpec.scala @@ -1,11 +1,10 @@ package symsim package examples.concrete.pumping -import laws.AgentLaws -import laws.EpisodicLaws +import symsim.laws.AgentLaws +import symsim.laws.EpisodicLaws class PumpIsAgentSpec extends SymSimSpec: checkAll ("concrete.pumping.Pump is an Agent", AgentLaws (Pump).laws) - checkAll ("concrete.pumping.Pump is Episodic", EpisodicLaws (Pump).laws) diff --git a/src/test/scala/symsim/examples/concrete/pumping/PumpSpec.scala b/src/test/scala/symsim/examples/concrete/pumping/PumpSpec.scala new file mode 100644 index 00000000..cd7853c5 --- /dev/null +++ b/src/test/scala/symsim/examples/concrete/pumping/PumpSpec.scala @@ -0,0 +1,71 @@ +package symsim +package examples.concrete.pumping + +import symsim.concrete.Randomized.given +import symsim.CanTestIn.given +import symsim.CanTestIn.* + +import org.scalacheck.{Arbitrary, Gen, Prop} +import org.scalacheck.Arbitrary.* +import org.scalacheck.Prop.forAll +import org.scalacheck.Prop.* +import examples.concrete.pumping.* + +// To eliminate the warning on Pumping, until scalacheck makes it open +import scala.language.adhocExtensions + + +class PumpSpec + extends org.scalacheck.Properties ("Pumping"): + + // Generators of test data + 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] (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) + + + // Tests + + property ("The tank level never can be negative") = { + 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 >= TankMin + } + } + + property ("There is an observable state for each valid state in the environment") = { + forAll (flow, head, head_mean, tank, time, water, past_head_mean, actions) { + (f, h, hm, tl, t, w, phm, a) => + val (s1, r) = Pump.step (PumpState (f, h, hm, tl, t, w, phm)) (a).head + exists (obsStates) {s => s == Pump.observe (s1)} + } + } + + property ("The time can be more than TankMax (comparing is not correct)") = { + exists (time) { t => + t > TankMax + } + } + + property ("Test getDemand") = { + forAll (time) { t1 => + exists (time) {t2 => Pump.demand (t1) != Pump.demand (t2)} + } + } + + property ("There is an action for each state which results no overflow") = { + forAll (flow, head, head_mean, tank, time, water, past_head_mean) { + (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 <= TankMax) + } + } + }