Compartmental System

This forum contains all archives from the SD Mailing list (go to http://www.systemdynamics.org/forum/ for more information). This is here as a read-only resource, please post any SD related questions to the SD Discussion forum.
Locked
"bjobj"
Junior Member
Posts: 4
Joined: Fri Mar 29, 2002 3:39 am

Compartmental System

Post by "bjobj" »

I should like to post preliminary responses to those that have answered to my plea for help with developing matrix models of compartmental systems.

Andy,

Thank you very much for your reply. Judging from the table of contents your new book seems to be very interesting and the appendix on spatial dynamics is certainly relevant. I shall try to get my hands on the book as soon as possible.

Colin,

Thank you very much for the Powersim equations for a compartmental matrix model. I would be thankful if you could send me the model so that I can see its structure.

I have not yet digested your comments on constraining the content of compartments with minimum and/or maximum values, but it does seem to be difficult to implement this feature in a Powersim model. Bob Eberlein has mentioned the possibility of doing this in Vensim. He says

"Now you need to figure out how to override the inflow
equation to prevent capacities from overflowing - that
requires allocation logic (Vensim has a function to do this
that might work)".

j-d,

You are right when you say that the problem is not well defined. I wanted to keep the model on a general level and use it as a kernel or module for more specific models. Your comments made me realise that it may be necessary to be more specific about the content of the compartments.

My point of departure was a system of compartments all of which contained various categories or kinds of the same "stuff". Let me take as an example a health service system, such as a hospital, which may be said to consist of a number of compartments, some of which may be physical while others may be conceptual. We may have patients waiting to be admitted, patients under treatment in various wards, patients waiting for discharge in various wards, patients that have been discharged, etc. Patients may "flow" through the system by being admitted, transferred between wards, transferred from the state of being under treatment to the state of waiting for discharge, and by being discharged. This could perhaps be an example of what you call case 1.1. It might be of importance to be able to assign maximum values to compartments in such a system, for example because a hospital ward has a finite size or a finite number of beds.

Each compartment in this system might contain different categories of patients, for example patients with schizophrenia, affective disorders, personality disorders, etc. The rate of flow through the system might be different for these different categories of patients. This could be an example of case 1.2.

I am in point of fact chief of just such a system - a fairly large acute psychiatric hospital with six wards.

Another example might be an epidemiologic system with people in various states of health or disease, and with various categories of people within each state. The rate of transfer between states might be different for different categories of people.

Thank your for the reference.

Kinds regards to all of you,

Oddur
From: "bjobj" <bjobj@online.no>
"bjobj"
Junior Member
Posts: 4
Joined: Fri Mar 29, 2002 3:39 am

Compartmental System

Post by "bjobj" »

Hi Geoff

Your points are well taken. I should like to discuss one of them. You ask why I am trying to build this vast model.

A model of a compartmental system does not have to be very large before a diagram of a scalar model becomes very complicated. Among the examples provided by Powersim with Constructor there are two models of a system of a population split into four age periods and two regions .One of the models is an array model called poparray.sim with the population being represented by an array in two dimensions. Births, deaths, and aging are modeled using flows into and out of the population array variable. I do not consider this a vast model. Nevertheles the scalar model called popsclr.sim of this system is very complicated. The design of the scalar model entails much more work than that of the array model and the diagram is much more complicated and even confusing. If you do not have access to these two models I can perhaps give you an idea of the difference between the models by comparing the diagram of the scalar model to a vast bowl of chinese noodles and that of the array model to an !
exquisitely simple japanese dish.

The diagram of the array model shows most of the important relationships between compartments, flows, auxiliaries and constants. It does not show the relationship between the individual compartments which is a drawback. This can however be compensated for by drawing the individual compartments and the flows between them with a program for drawing flow diagrams. It is much easier to design a user interface for an array model of a comparmental system than for a scalar model and user interaction is also easier.

Assigning minimum and/or maximum values to the content of compartments in a scalar model will make even a model of a relatively small compartmental system very complicated. I therefore hope that it will prove possible to design a matrix model of such a system.

Regards,

Oddur
From: "bjobj" <bjobj@online.no>
"bjobj"
Junior Member
Posts: 4
Joined: Fri Mar 29, 2002 3:39 am

Compartmental System

Post by "bjobj" »

I have used Powersim Constructor 2.51 to implement an array model of a compartmental system with flows possible between all compartments. I am sure somebody has designed a more efficient, more elegant, or in other ways better model. Can somebody point me towards work that has been done in this area?

I have also tried to extend the model so that a maximum and minimum value can be assigned to each compartment - so far without success. Is anybody aware of such models?

[Hosts note: I asked Oddur to elaborate a little]

Matrix model of a system of n compartments with intercomparmental flows.

In order to make my questions more specific I take as my point of departure a system with four compartments.

Let the 4 compartments of the system be denoted by C_1, C_2, C_3, C_4. Let the content of compartment C_i be denoted by c_i. Let the matrix C = (c_i).

Let the rate of flow per unit time from C_i to C_j be denoted by f_i,,j. Let the matrix F = (f_i,,j).

Let the rate of flow per unit content per unit time from compartment C_i to compartment C_j be r_i,j. Let the matrix R = (r_i,j).

Let f_i,j = r_i,j*c_i.

How can a system of this description be represented by a matrix model in Vensim, Powersim or Stella Research?

How can the system be represented by a matrix model if it has to allow for the possibility of assigning minimum and maximum values to the content of the compartments?

How can flow from a source to the compartments and from the compartments to a sink be incorporated into the matrix model?
Colin Brown
Junior Member
Posts: 2
Joined: Fri Mar 29, 2002 3:39 am

Compartmental System

Post by Colin Brown »

Oddurr,

I tried to do something like this several years ago, modeling the flow
of customers between different competing companies. It is relatively
easy to represent the basic system you described using SD modeling
software. Adding additional sources and sinks is fairly straightforward
too. This problem gets difficult when you attempt to add min and max
constraints onto your containers. I dont believe there is a way to do
this without solving a set of simultaneous equations (creating circular
references...) to handle the flows when a constraint is encountered
(though if anyone does know how to do this, please let me know!).

For example, using your notation, we have two containers C_1 and C_2,
the flow from C_1 to C_2 is F(1,2) and the flow from C_2 to C_1 is
F(2,1).

If C_1 - F(1,2) + F(2,1) > Max_Value then C_1 is constrained.

We could then say that G(2,1) is the flow from C_2 to C_1 which is not
constrained, and H(2,1) is the flow from C_2 to C_1 which is
constrained, so F(2,1) = G(2,1) + H(2,1), and G(2,1) = Max_Value - C_1 -
F(1,2). So F(2,1) is really the "desired flow", G(2,1) is the "actual
flow".

So we have solved the problem of C_1 possibly going over its maximum
value, but have created H(2,1), the amount of the original flow from C_2
to C_1 that "could not flow" due to capacity constraints. In order to
keep material from creating/destroying itself in the system, you now
need to go back and adjust F(2,1), because only G(2,1) is really
flowing. So if we try to set the outflow F(2,1) to G(2,1) we are
creating a circular reference. You really need to do some kind of
iterative algorithm to converge on the true "actual flow" each timestep.

So, G(2,1) = Max_Value - C_1 - F(1,2) = F(2,1) - H(2,1), and vice-versa
for G(1,2), giving us the set of simultaneous equations:

H(2,1) = F(2,1) - Max_Value + C_1 + F(1,2)
H(1,2) = F(1,2) - Max_Value + C_2 - F(2,1)

The numerical problem becomes starting with the "desired" values of
F(1,2) and F(2,1) (computed from multiplying the matrix R by the
container vector), and iterating them to find the values that make
H(2,1) and H(1,2) equal zero. Check out any numerical analysis textbook
for more detail on this subject.

If youre using Powersim you could write a VB or C++ function that sits
on top of the model and does this, though if this is all your model
needs to do, you are probably better off just writing a computer program
to do the entire simulation.

There are ways you could make your model solve these simultaneous
equations through simulation, but they are somewhat ugly. If you did it
all in one model, you would have to make your timestep incredibly small.
This is because the simultaneous equation algorithm may need 100s of
timesteps (or more!) to converge within a specified error tolerance, and
the rest of your model will be simulating away the whole time. This
also creates the need to have an array element for each timestep in the
simulation to store intermediate iteration values because a new set of
equations to solve will be created every timestep. Very Ugly!

If you use Powersims Co-model feature to do the iterations in a
separate co-model running on a different timescale (something like 1000
timesteps for each timestep in the main model) you might be able to do
this. This would eliminate the problem of having multiple iterations
going on at the same time and the need for an array element for each
timestep. The only drawback I see with this is that there is a 1
timestep delay for data transfers between co-models. This would only
generate negligible error if your main model timestep is sufficiently
small.

In summary, unless other aspects of your simulation make you really want
to use system dynamics software for this model, you might want to
consider just writing a computer program, or writing a program in some
math packages programming language like Matlab.

If you can live without the min/max constraints, here are the Powersim
equations to do this elegantly. I can send you this model too if you
like.

dim Compartments = (i=1..4)
init Compartments = [5, 10, 2, 4]
flow Compartments = -dt*Sink_Outflow
+dt*Source_Inflow
+dt*Inflow_j
-dt*Outflow_i
dim Inflow_j = (j=1..4)
aux Inflow_j = ARRSUM(F(1..4,j))
dim Outflow_i = (i=1..4)
aux Outflow_i = ARRSUM(F(i,1..4))
dim Sink_Outflow = (i=1..4)
aux Sink_Outflow = .1*Compartments
dim F = (i=1..4, j=1..4)
aux F = R*Compartments(i)
dim Source_Inflow = (i=1..4)
const Source_Inflow = 1
dim R = (i=1..4, j=1..4)
const R = [[0, .25, .25, .25],
[.33, 0, .33, .33],
[.1, .1, 0, .1],
[.2, .2, .2, 0]]

I made a constant "source_inflow" of 1, and a "sink_outflow" of .1 times
the compartments current value, chose some initial values and made up
some fractions for the matrix R.

Hope this is helpful,

Colin Brown
cbrown@carlisle-co.com
978-318-0500 ext. 179
"j-d"
Junior Member
Posts: 5
Joined: Fri Mar 29, 2002 3:39 am

Compartmental System

Post by "j-d" »

This is a very interesting problem. I have a few comments and then some
suggestions as to how you can solve it, with a reference.

(1) First, it seems to me that some elements of the problem are still not
well-defined. Is there only one product flowing from one compartment to the
other or are there >1 products, i.e., does each compartment send out a
different product? To illustrate from your example of 4 compartments:

Case 1.1: A, B, C, and D are 4 "compartments" (possibly stocks) which all
carry apples and apples flow out of A, B, C, and D at the rate dependent on
the respective stocks of apples in each of A, B, C, or D;

Case 1.2: A, B, C, D are 4 "compartments" (possibly stocks) which
respectively carry positive "A"pples, "B"ananas, "C"herries, and "D"ates, at
time t=0, and these four flow out of A, B, C, and D respectively at the
rates dependent on stocks of "A"pples, "B"ananas, "C"herries, and "D"ates,
respectively. At time t>0, A, B, C, D each will likely contain all four
products!!!

>From your problem definition, it seems A could send out stuff to A, B, C, or
D, for each of the above cases, and so on for B, C, and D. In case 1.1, the
problem, at least to me, is not very interesting because after a sufficient
time, you will have some apples (+ve, -ve, or 0) left in each stock, but
because all apples are identical, you wont know which apple came from which
stock, and also the calculations and model are pretty straightforward there.
So I will assume you had Case 1.2 in mind, which is very interesting and
complicated. You now have 4 ARRAY stocks, A, B, C, and D, each of dimension
4 ("A"pples, "B"ananas, "C"herries, and "D"ates) and each are sending stuff
to each other. You can start with

A = (100, 0, 0, 0) --> 100 apples, ZEROs rest
B = (0, 200, 0, 0) --> 200 bananas, ZEROs rest, etc., etc.

Now the neat thing happens: at the start of the simulation, A "pushes" its
apples to B,C, D, and B pushes its bananas and so on - after some time, all
of A, B, C, D have some (+ve, -ve or zero) amount of each product and we
dont know how the system will behave after this, because their rates of
"push" depend on the quantities in stock of the "pusher" products and now
ALL dimensions in ALL stocks are active. I havent tried this in Powersim,
but it may not be easy.

There is a reason why "push" is in quotes. I dont know the context of your
problem, but I was intrigued that the flow f_i,j depended only on the
content c_i of C_i. In real life, one would expect flows to depend on some
kind of potential difference, pressure difference, height difference,
voltage difference, supply and demand mismatch, etc, etc. Your case reminded
me of the "push technology" in todays networked economy - "Whether you want
this news or ad or not, we will "push" it into your e-mailbox". So - what
problem are you working on? A context could also tell us how to go about its
solution.

Also, is there a purpose for this implementation? I mean, minimizing costs,
flow times, etc. or are you just simulating some observed network behavior?

Most important point - the above problem is a non-typical system dynamics
problem. Since you mentioned "n" compartments and "matrices", I would say
why dont you go the route of operations research, where problems much worse
than this have already been solved using matrix algebra. Implementation in
system dynamics is fine and fun, but if you want to solve the problem
quickly, here is the reference:

"Algorithms for Network Programming" by Jeff L. Kennington and Richard V.
Helgason, John Wiley & Sons, 1980 edition - chapter 6, pages 124-165. This
chapter deals with "The Multicommodity Network Flow Problem" as I observed
in Case 1.2 above. In this case it is an optimization problem over the
network, for example, there would be demands for ("A"pples, "B"ananas,
"C"herries, and "D"ates) at each node, corresponding prices, different rates
of product transportation costs and so on, and we want to minimize the total
costs. Also it will be a static solution - at one point in time - to make it
dynamic optimization, I would suggest writing a separate program (use
dynamic programming for this, though let us not get into that now!!). Loads
of your best Scandinavian coffee will help too.

Hope I havent made the problem much more complicated than you had in mind.
If you need specific examples/solutions, let me know and I can put them on
my webpage (doing math in a public forum using ASCII is criminal).

Jaideep

Jaideep Mukherjee, Ph. D.
From: "j-d" <j-d@technologist.com>
713 523 2713; www.netopia.geocities.com/shunya/
"Phil Odence"
Junior Member
Posts: 15
Joined: Fri Mar 29, 2002 3:39 am

Compartmental System

Post by "Phil Odence" »

[the model is located at http://www.vensim.com/sdmail/models/SD2339.ITM ]

The attached ithink model (can be opened with STELLA; you can get a free
runtime at www.hps-inc.com, if you dont have the software) handles Oddurs
compartment question rather compactly (4 elements). The trick is to
calculate how much material should flow from and to the one dimensionally
arrayed stock in a two dimensional matrix, and then to sum the outflows
(rows) and the inflows (columns) and flow the difference.

How to implement mins/maxs depends on particulars of how the real system
works. When it a compartment gets full does it shut off inflows, increase
outflows, both? Even once this is specified, it can get tricky to work out
the logic, and you have to be careful to conserve material in that I did not
model the flows as conserved.

External inflows/outflows is quite straight forward.

Phil Odence
High Performance Systems
From: "Phil Odence" <podence@hps-inc.com>
kevin.a.agatstein@us.arthurander
Junior Member
Posts: 6
Joined: Fri Mar 29, 2002 3:39 am

Compartmental System

Post by kevin.a.agatstein@us.arthurander »

Oddurr,

In regards to the solution put forward by Colin Brown in response to his
question, there may be a formulation error in the following equation and wish
discussion:

If C_1 - F(1,2) + F(2,1) > Max_Value then C_1is constrained
where Cs are containers (or stocks)
and Fs are flows

By definition, you cannot substract a flow from its stock, as for one there is a
dimensional inconsistancy. The flows into and out of a stock have the unit of
measure of the stock over time. Hence, you cannot substact units per time from
units (people per month from people for example.) Also, because of this, I
question whether Colins solution is robust over any reasonable solution
interval (TIMESTEP or DT), and not just a solution interval of 1.

I do believe though for the most part Colins proposed solution to be quite
interesting though.

One question you need to address when considering the "constraint" is what to do
when a compartment is below its constraint but the next step will kick it over,
do you want the compartment to fill to its constraint or to not fill at all.
For the former, you need to pulse in the difference of what will be in the stock
minus the maxium. In the later case, simply "turn off" the inflow at the
constraint.

Best of luck with your modeling efforts,
Kevin Agatstein
Arthur Andersen, LLP
From: kevin.a.agatstein@us.arthurandersen.com
Bill Braun
Senior Member
Posts: 73
Joined: Fri Mar 29, 2002 3:39 am

Compartmental System

Post by Bill Braun »

Phil wrote:

>An example (that is eluded to in the ithink Business Applications Guide):
>When hospital administrators reviewed a map theyd created with nurses in
>the emergency care facility, they discovered that unbeknownest to them there
>was as significant flow of patients being turned away on a daily basis due
>to an overcrowed waiting room. The issue would never have come up in the
>discussion had the stock of patients in the waiting room been buried as
>"C_3" in an array.

I find myself in precisely this dilemma regards to a model I am working on.
The complexity of the model makes the use of dimensions and arrays very
attractive - it certainly simplifies the modeling process itself. As Phil
points out, it also makes it near impossible for the non SD person to make
any sense of it and gain any insights and learning.

Currently I am working along the line of buidling a GUI in Paradox (through
a DDE link) that will "translate" the model and make the equivalents of the
"C_3" stock (which Phil made on example of above) visible and meaningful to
the user. Good example of the prinicple of tradeoffs: I gained some ease in
modeling and gave up some ease in its value as a learning tool.

Bill Braun
From: Bill Braun <medprac@hlthsys.com>
Locked