Discrete Control tests based on realistic model

[SouthEndMusic]

At the time of writing controlling node types other than Pump has been merged (#424), and making nodes active and inactive is on its way (#425). These features do not have proper tests yet, but they came about form a discussion about a model of the Netherlands. So, this issue is for

• Constructing said model;
• Testing discrete control on non-pump nodes and node activity.
• Tests for PID controlled pump #482

For the details of this model see below.

[visr]

[gijsber]

FlowBoundary holds a timeseries ranging from 200-1000m3/s with some variation around the 250 and 800m3/s
DiscreteController 17 and 18 activate node 8 and node 13 and deactivate node 9 and 14 when the FlowBoundary > 800m3/s. This state is preserved till the FlowBoundary drops below 750 m3/s.

If Flow Boundary between 250 and 800m3/s, the flow rate of pump 9 is set to 25m3/s by controller 18 (using DC or PID). If Flow Boundary below 250m3/s, the flow rate of pump 9 is set to 15m3/s. This flow rate is held until the Flow Boundary rises above 275m3/s.
PID control 20 preserves the water level of basin 12 at 6m+NAP.

Basin 2 (length 10km) has the following characteristics. All other basins are in the same ballpark

level	totalwidth
1.86	0
3.21	88
4.91	137
6.61	139
8.31	141
10.07	219
10.17	220
10.27	221
11.61	302
12.94	606
13.05	837
13.69	902
14.32	989
14.96	1008
15.59	1011

[SouthEndMusic]

@gijsber thanks for the information. I have some questions:

• Are the areas equal to the totalwidth times the length?
• What are the parameters of LinearResistance (3, 4, 11), TabulatedRatingCurve (8, 13), PidControl (20), LevelBoundary (16), Pump (14)?

Other points:

• Time dependent flow boundaries are not yet supported: Time-dependent flow boundary #439
• See here for the current draft script of the model: Ribasim/python/ribasim_testmodels/ribasim_testmodels/dutch_waterways.py (branch Dutch_waterways_model)

[gijsber]

Are the areas equal to the totalwidth times the length?

Yes.
Here are some best guesses (actually this is a task for Hydrologic):
LinearResistance: 5000 (assumption)
LevelBoundary: 3 (assumption)
Pump: min 0, max 50, default setpoint 25 (assumption)
PID: Kp: -0.005, Ki:0, Kd: -0.002 (from SBK3 model for gate)
QH(8): 7.45m Q=418m3/s
elke centimeter stijgt Q met 2.15m3/s
QH(13) is vergelijkbaar maar dan 3 meter lager

[SouthEndMusic]

@gijsber I have distilled the discrete control of nodes 17 and 18 as follows.

Notation

• q: Flow at edge (1,2), flow of flow boundary with id 1
• w: Flow at edge (6,8), used to indicate whether node 8 is active

Conditions

1. q > 250
2. q > 750
3. q > 800
4. w > 0

Control states

• "pump_low":
  • pump (node 9) active
  • pump (node 9) flow rate 15
  • TabulatedRatingCurve (node 8) inactive
• "pump_high":
  • pump (node 9) active
  • pump (node 9) flow rate 25
  • TabulatedRatingCurve (node 8) inactive
• "rating_curve":
  • pump (node 9) inactive
  • TabulatedRatingCurve (node 8) active

Logic
To massively reduce the number of logic table entries needed, I want to introduce the wildcard A for any as an addition to T and F for true and false respectively.

• FAAA -> "pump_low"
• TFAA -> "pump_high"
• ATFF -> "pump_high"
• ATFT -> "rating_curve"
• AATA -> "rating_curve"

Final notes

• I am aware that using the flow condition on edge (6,8) is an ugly thing to do, but the problem is that this 'memory for control' that we discussed for 'being on your way to a setpoint' is a different type of condition than the ones currently supported and fitting in the framework of VectorContinuousCallback.
• This requires an introduction of the wildcard A: Add truth state wildcard for DiscreteControl #440.
  Also possibly some validation on the logic data (no effective truth state duplicates). Note though that not all truth states have to be covered because many are impossible, for instance FTAA.

[gijsber]

Hi @SouthEndMusic, That looks pretty good. We miss one condition (q>275). If we add that it would become:

Conditions

1. q > 250
2. q > 275
3. q > 750
4. q > 800
5. w > 0

Logic

FAAAA -> "pump_low"
TFAAT -> "pump_low"
TFAAF -> "pump_high"
AATFF -> "pump_high"
AATFT -> "rating_curve"
AAATA -> "rating_curve"

[SouthEndMusic]

I realized a potential problem with listening to a flow influenced by a time series as a condition, and that is that such flows are not continuous. This might cause problems in the root finder of VectorContinuousCallback. Thoughts @visr?

[visr]

Yeah I think indeed that would be a problem, though I guess we could test to see what it does in that case. If there are errors we should probably go for a DiscreteCallback.

[SouthEndMusic]

What about the option of linear interpolation of the time series? Though that makes integration a bit more complex.

[visr]

Yeah I've gone back and forth on that, see #188 for the reasoning for the current approach. Though the discontinuities in the flow are from the callback, and I think that is allowed to do that there, so perhaps the ContinuousCallback can handle it fine, if it only assumes continuous functions between tstops?

I've looked a bit at similar arguments online but couldn't find much more than SciML/ModelingToolkitStandardLibrary.jl#123 (comment)

[gijsber]

From a hydrological perspective, I would think that linear interpolation is acceptable.

[visr]

I think it depends, e.g. if you have a non-equidistant precipitation timeseries like

2020-01-01 00:00:00 0.0
2020-02-01 03:02:40 10.0
2020-02-01 03:06:00 0.0

which represents a short intense shower for a few minutes, do you want to gradually go from 0 to 10 over the month of January?

If it is equidistant it might be good enough in general.

[SouthEndMusic]

@gijsber in the image you sent of the model edge (5,7) to the terminal is invalid, there must be a flow prescribing node in between.

[gijsber]

For precipitation timeseries, this assumption is indeed not wise, but I would think that flow series have more continuous behaviour (and indeed you provide them typically as equidistant series)

[gijsber]

@SouthEndMusic I would suggest to add a linear resistance as needed

[SouthEndMusic]

@SouthEndMusic I would suggest to add a linear resistance as needed

A LinearResistance is not allowed between a Basin and a Terminal because the Terminal has no level. Shall I replace the Terminal with a LevelBoundary with an assumed level?

[gijsber]

using a LevelBoundary is fine to me

[SouthEndMusic]

Using #447 and #440,

With conditions

1. q > 250
2. q > 275
3. q > 750
4. q > 800,

I think this logic would work:

• FAAA -> pump low
• UAAA -> pump low
• TAFA -> pump high
• ATAF -> pump high
• AAAD -> rating curve
• AAAT -> rating curve.

The 'pump high' cases have some overlap which should be cleaned up.

[gijsber]

Suggested initial levels:
basin 2: 8.31m+NAP
basin 6, basin 5: 7.5m+NAP
basin 10: 7.0m+NAP
basin 12: 6m+NAP
basin 15: 5.5m+NAP

[SouthEndMusic]

@gijsber with the above initial levels the simulation now converges but only without the analytical Jacobian. The result doesn't seem very sensible though:

[SouthEndMusic]

With much lower resistances (1e-2 instead of 5000) the result already looks better:

[SouthEndMusic]

This at least explains why the PID controlled basin doesn't reach its target level: the PID controlled pump can only take out water while the level is generally too low: