Enforced (Stella Style?) Limits and Warning (Vensim Style) Limits

[gherka]

When I set limits to stocks (say, 0-1000), PySD doesn't seem to respect them and I get runaway values in my SIR model. Is there some other setting I'm missing? When I run the model in Stella it works correctly. model.doc() shows that limits are set for stocks, but that doesn't seem to have any impact. I've attached the .xmile model from Stella.
SIR model.zip

[JamesPHoughton]

Hi @gherka,

There are two schools of thought as to what a limit should do in a model - one school says that a limit should be a hard boundary, and if a stock was going to exceed that limit, then the modeling software should force no more growth of that stock. I believe this is what Stella does, and is probably why your model works in that environment. The other school says that the limits shouldn't modify the behavior of the model, but should throw a warning to the user to let them know that the stock has exceeded its limit. This is what Vensim does.

The logic of the second school is that the structure of the system itself should be explicit about constraining the value of the stock. Here's a toy example - say that a room has a capacity of 100 people. At first the inflow rate is constrained by the number who can get through the door, and so the population in the room grows linearly. Then as the room starts to become full, the rate limiting factor is how fast people can navigate the room and find seats. So, you get gradual adjustment to the capacity. Here's what that would look like over time:

Stella lets you just run up against the stop at full speed, which is fine in a lot of cases because it won't matter to the question your model is trying to answer. Vensim says you have to be explicit about the behavior you want as you get close to the constraint, and if you want a hard stop vs a soft landing you have to build that behavior into your model. Here's an example for a soft landing:

Characteristic Time for Occupancy Constraint=
		5
	Units: Minute
	
Max Occupancy=
	100
Units: Persons

Maximum Inflow Rate Given Current Occupancy=
	(Max Occupancy - People in Room)/Characteristic Time for Occupancy Constraint
Units: Persons/Minute [0,?]

Nominal Inflow Rate=
	5
Units: Persons/Minute

People Entering=
	MIN(Maximum Inflow Rate Given Current Occupancy, Nominal Inflow Rate)
Units: Persons/Minute

People in Room= INTEG (
	People Entering,
		0)
Units: Persons

In PySD, I've chosen to take the Vensim school, mostly because I started with Vensim and know it best. We might in the future add a hard_limits option, but for the time being, you'll probably be best off writing the constraint into your model.

I've been working on some bounds testing functions that are post-simulation tests of behavior. They aren't documented yet, but they are here in the code. I'm afraid something seems to be broken at the moment (in the current repo) but eventually you should be able to do:

model = pysd.read_xmile('SIR model.xmile')
model.doc()

Real Name	Py Name	Unit	Lims	Type	Eqn	Comment
Dead	dead	People	(0.0, 1000.0)	component	daily_deaths
Infectious	infectious	People	(0.0, 1000.0)	component	sick_per_day - daily_recoveries - daily_deaths
Recovered	recovered	People	(0.0, 1000.0)	component	daily_recoveries
Susceptible	susceptible	People	(0.0, 999.0)	component	- sick_per_day
daily deaths	daily_deaths	People/Days	(None, None)	component	Infectious*(1-survival_rate)/disease_duration
daily recoveries	daily_recoveries	People/Days	(None, None)	component	Infectious*survival_rate/disease_duration
disease duration	disease_duration		(None, None)	constant	5
infection rate	infection_rate		(None, None)	constant	.3
sick per day	sick_per_day	People/Days	(None, None)	component	InfectiousSusceptibleinfection_rate
survival rate	survival_rate		(None, None)	constant	.25

and then:

mat = pysd.testing.create_bounds_test_matrix(model)
mat

	Real Name	Comment	Unit	Min	Max
Dead		People	-inf	inf
Infectious		People	-inf	inf
Recovered		People	-inf	inf
Susceptible		People	-inf	inf
daily deaths		People/Days	-inf	inf
daily recoveries		People/Days	-inf	inf
disease duration			-inf	inf
infection rate			-inf	inf
sick per day		People/Days	-inf	inf
survival rate			-inf	inf

(There are clearly some errors in here at the moment!)

and then:

pysd.testing.bounds_test(model.run(), bounds=mat)

unfortunately when I do this at the moment, I'm getting NaN, so there is some debugging to do.

Sorry not to be more help - hope you work out what you need.

James

[JamesPHoughton]

I'm going to leave this open as a feature request for future dev, so that we can put in an option for hard limits down the line. @alexprey, you probably know more about this than I do. What does the XMILE spec have to say about the different types of limits?

[gherka]

Thanks a lot for your explanation, James. Very helpful. I wonder if the issue in my example is not just about limits, but also non-negative stocks. When I created the stocks in my model, I set them to be non-negative and that's recorded in the .xmile file. However, PySD coversion is happy to have "negative" number of people in the Susceptible stock under certain outflow conditions. I notice that in your equations you explicitly defined what should happen in case of a negative stock:

People in Room= INTEG (
	People Entering,
		0)

I couldn't find an equivalent function in Stella, except for setting the stock as non-negative.

Perhaps the code generated when you load the model could check for that? Assuming it's possible to pass the non-negative flag from .xmile parser downstream...

@cache('step')
def susceptible():
    """
    Real Name: Susceptible
    Original Eqn:  - sick_per_day
    Units: People
    Limits: (0.0, 999.0)
    Type: component
    Non-negative: True #new comment line
    """
    if non_negative & (integ_susceptible() < 0):
        return 0
   
    return integ_susceptible()

integ_susceptible = functions.Integ(lambda: -sick_per_day(), lambda: 999)

[JamesPHoughton]

Ah, no you've given me too much credit. In the room capacity model, I make the assumption that there is no outflow. (The 0 in the INTEG equation sets the initial condition for the stock). If I want to enforce that there be no negative level in the stock I would probably make something like this:

Average Time in Room=
		100
	Units: Minute [0,?]
	
People in Room= INTEG (
	People Entering-People Leaving,
		0)
Units: Persons

People Leaving=
	People in Room/Average Time in Room
Units: Persons/Minute [0,?]

The first-order control on the outflow is essential to make sure the flow doesn't go negative. If I thought the room exit rate had a limited flow rate, I could do a formulation like I have on the inflow, with the MIN function switching between whichever condition was limiting, the number of people who want to leave, or the maximum outflow rate.

With the complete model I end up with the inflow and outflow coming to equilibrium somewhere below the max capacity:

[alexprey]

Hi every one! 😸
@JamesPHoughton by xmile we have range node for stocks and flows, but this configuration has effects only for input values. But we have in stocks - non_negative value. And this behaviour can be defined for whole model. Also, xmile have interesting options to change it behaviour - queue, conveyors, leaking.

But, xmile is extendable language and I think it's a good point to think about introducing additional attributes and tags for stocks and flows to define custom range of value bounds.

---

Regards, Alexey Mulyukin
Lead core-developer of sdCloud.io dev team.