Optimal Dynamical Balance of Raw Materials —
Some Concept of Embedding Optimization in Simulation
on System Dynamics Models and Vice Versa
Elzbieta Kasperska, Elwira Mateja-Losa, Damian Stota
Institute of Mathematics
lesian University of Technology
ka 23, Gliwice 44-100, Poland
e-mail: elakaspe@polsl.gliwice.pl, elwimat@polsl.gliwice.pl, damslota@polsl.gliwice.pl
Abstract
The purpose of this paper is to present the optimal dynamic balance of raw materials in
two formulations. The model DYNBALANCE(3-1) appears in two cases, named I and II.
The main idea of case I of model DYNBALANCE(3-1)concist in optimal control of sup-
plying three raw materials according to productivity of three technologies, by taking into
consideration possibilities of sources of raw materials and actual demand for product
on market. This optimal control lies in optimal balancing of raw materials, according to
plans: cost of raw materials, cost of production and according to forrecast of demand. The
obj function in optimizing process (in Coyle’s sense) is consisted of three main com-
jecti
ponents and one component which has ” penalty” function. The optimization experiments
were performed usying COSMOS package (Computer — Oriented System Modelling Opti-
misation Software), which is software tool automatically linking a dynamic simulation
model to an optimization package. In this paper authors present the whole structure of
model DYNBALANCE(3-1) case I, with the legend and list of variables and parameters.
The equations of the model reader can find in program output in Appendix B.
The main idea of case IL of model DYNBALANCE(3-1) concist in replecement of
problem of solving of balance of raw materials, which minimize the norm ||Ax — b|| (at
the condition x > 0) by the problem of solving pseudosolution of extended balance of raw
materials, which minimize the norm |{Ax — || (without condition x > 0). The benefit
of such solution of problem is that we can embedding optimization in simulation on
System Dynamics model by using language Dynamo. In each step during simulation the
program solves the x by the Legras formula: (Z’A)~!A"S, step by step. Such a solution
(balance of raw materials) we can named “dynamic”, so the embedding the optimization
in simulation, on System Dynamics models, has dynamical character.
Many interesting experiments were undertaken by authors and discussed in paper. At
the end of paper some conclusions were formulated, specially concerning the problem of
simplificationes of models.
Keywords: System Dynamics Method, Optimization Embedding in Simulation, Si-
mulation during Optimization, Balance of Raw Materials.
1 Introduction
The idea of extending System Dynamics method was undertaking already by authors in pa-
per [5]. The point was that some known methods of optimization can be embedding in simu-
lation on System Dynamics model and vice versa. The authors were occupied with dynamical
balance of production. Now, in this paper, the problem is opposite. The balance relates to
three raw materials and one product. The two cases of model, called DYNBALANCE(3-1),
were created by Kasperska and many interesting experiments were performed using Profes-
sional Dynamo 4.0 and COSMIC and COSMOS by Mateja-Losa and Stlota.
the sum function of
fitting demand (sffd)
the balance
of raw
materials
level of inventory level of a material during
of production (lin) oO transformation (Imt)
Source 2
raw material
< Q tchni3
the sum function
of fitting cost balance
of production (sffpr) i ovum
the sum function of fitting
cost balance of raw
materials (sffrm)
the sum loose
of profit (lopr)
Figure 1. The optimal dynamics balance of raw materials (model DYNBALANCE(3-1))
2 Some examples of the optimal balance of raw materials
In this paper we present two cases of model DYNBALANCE(3-1). Generally speaking, in
the case I the model has the structure, which allow optimization in sense of Coyle [2-4]. In
the case II the model has the structure which allow embedding optimization in simulation on
System Dynamics model, in the sense of authors.
Lets, at first, present on Figure 1 the structure of model DYNBALANCE(3-1) in case I.
The graphic convention is based on Lukaszewicz’ symbols [7,8], supplemented by own idea
of authors, and is explained in legend to the figure.
Legend to the Figure 1
Entering flow
sical level in Forrester sense;
Level
Leaving flow
ical rate (action) in Forrester sense;
Rate
ical source (independent variable);
flow of information;
flow of materials;
A
om 2 direct amplification (A — amplification coefi-
cient);
R=AlL
symbol in double frame (local relationship in
model, for example the optimal balance of raw
material);
oe
rR parameter;
Q
, N
AY resultat
LU OR vaviable
auxiliary variable.
—. mathematical ¥
pa operation
oe Ts.
The main idea in case I of model DYNBALANCE(3-1) consist in optimal control of
supplying three raw materials according to productivity of three technologies, by taking into
consideration possibility of source of raw materials and actual demand for product on market.
Optimal control lies in optimal balancing of raw materials, according to plans: cost of raw
materials, cost of production and according to forrecast of demand. The objective function
‘isted of three main components and one
in optimizing process (in Coyle’s sense [2-4]) is coi
component which has ” penalty” function. Lets explain the variables and parameters, to clear
the role of objective function.
Variables and parameters of model DYNBALANCE(3-1)
Levels:
Imt
lin
sffd
sffrm
sffpr
level of material during transformation;
level of inventory of production;
level which sum up the discrepancies between the mass plan of pro-
duction (equal to forecasted DEMAND ~ independent variable in that
model) and its realization (production are product from three raw ma-
terials, according to three technologies) , this level authors named the
sum function of fitting demand”;
level which sum up the discrepancies between the planed cost of raw
materials (total cost plan: parameter term) and its realization (which
depends on unitary cost of raw materials and actual calculated values
of three flows of raw materials);
level which sum up the discrepancies between the planed cost of pro-
duction (technology and peaples) (total cost plan: parameter tepr) and
its realization (which depends on unitary cost of production from three
technologies and actual calculated values of three flows of raw materials);
level which sum up the loose of profits, when actual sale is lower then
demand.
rate of raw material (in unit of production) which enter the level of
material during transformation, is calculated by appropriate summing
of optimised three flows of raw materials;
rate of raw material number 1 (which depends on optimal productivity
of first technology and actual cindition of source of raw material 1);
rate of raw material number 2 (which depends on optimal productivity
of second technology and actual cindition of source of raw material 2);
rate of raw material number 3 (which depends on optimal productivity
of third technology and actual cindition of source of raw material 3);
rate of production (we assume that level of production (material during
transformation) is delay of first order);
rate of sale (depends on actual rate of demand and possibility of inven-
tory of production);
rate of discrepancies between the mass plan of production and its actual
realization (input to level sffd);
rate of discrepancies between the cost plan of raw materials and its
actual realization (input to level sffrm);
rate of discrepancies between the cost plan of production and its actual
realization (input to level sffpr);
rate of loose of profits, when actual sale is lower then demand (we assume
the constant price of product).
Auxiliaries:
sourcel
source2
source3
rd
frd
dod
fob
Parameters:
ql
q2
g
tchn1
tchn2
tchn3
tpr
po
pl
perd
wl
w2
ws
p2
ucr1
ucr2
ucr3
ucpr1
ucpr2
ucpr3
term
tepr
cena
kara
actual possibility of supply of source of raw material 1;
actual possibility of supply of source of raw material 2;
actual possibility of supply of source of raw material 3;
demand for product (we assume sinusoidal character of curve of de-
mand);
forecasted demand for product (we assume sinusoidal character of curve
of forecast);
difference of demand from sale (measures the difference between actual
demand and its realization (sale));
objective function (sum of sffd, sffrm, sffpr and penalty which me-
asures the loos of profits, when the actual sale is lower then demand;
The components of the sum have weight factors.
product multiplier from the first raw material;
product multiplier from the fsecond raw material;
product multiplier from the third raw material;
ivity of first technology;
ity of second technology;
productivity of third technology;
time of production;
parameter step of input (demand);
parameter of amplitude of sinusoidal input (demand);
parameter period of input (demand);
weight factor for component sffd of function fob;
weight factor for component sffrm of function fob;
weight factor for component sffpr of function fob;
parameter of amplitude of sinusoidal characteristic of forecast of demand
(frd);
unitary cost of raw material 1;
unitary cost of raw material 2;
unitary cost of raw material 3;
unitary cost of production from raw material 1 (technology 1);
unitary cost of production from raw material 2 (technology 2);
unitary cost of production from raw material 3 (technology 3);
total cost of raw materials (plan);
total cost of production (technology and people) (plan);
price of product on market;
coeficient of gain of loose of profit (weight factor for lopr) (see: objective
function).
The optimization experiments were performed usying COSMOS [1]. COSMOS (the Com-
puter — Oriented System Modelling Optimisation Software) is a software tool which automati-
cally links a dynamics simulation model to an optimisation package. The type of optimisation
we used were so called ” Direct Optimisation”, which requires:
an objective function,
a group of parameters to be searched,
the permissible range for each parameters.
The objective function ia a variable which summarise the undesirable characteristic of
the system. The parameters are usually constant in the model. The author of COSMOS
package, prof. G. Coyle, didn’t discover which optimisation technique were used in package.
We suppose that it could be, for example, method COMPLEX for nonlinear optimisation. To
apply COSMOS we not need the knowledge of method but only some rules of optimisation
dialoque, which has form like on Figure 2.
Name of Objective Function
Minimise or maximise?
Parameters and ranges
Optimisation control
Parameter value printed
Length of simulation
Number of iteration
Step multiplier
Number of output lines
Type of optimisation
Base vector
Simplifier
Planning horizon
If none of these them
Direct optimisation
Main command mode
or
Rerun mode
Figure 2. General form of optimization dialoque (from Coyle [1])
Now, lets pay attention on case II of model DYNBALANCE(3-1). Some of the variables
and parameters don’t require the explanation (see: list of variables and parameters of model,
case I), but there are some new elements (see Figure 3). Theirs roles are easily visible on
Figure 2 and in Appendix B.
level of inventory
of production (lin)
level of a material during
transformation (Imt)
ucprl,2,3
b4,5,6,
PAP)
uerl,2,3
suml
level sum up the square
of discrepancies of first
equation (Ax—b)
level sum up the square
of discrepancies of third CS
equation (Ax—b)
sum2 summ
level sum up the square of] aoa
C) discrepancies of second @----~
equation (Ax—b)
Figure 3. Structure of model DYNBALANCE(3-1) ~ case II (with the elements of optimization in sense of authors)
The main idea of the method used is the replacement of the one problem named (1) by
the second problem named (2). The problem (1) consists in minimizing the norm ||Az — 6]|
by the condition « > 0. The matrix A has dimensions 3 x 3, the matrixes x and b have
dimensions 3 x 1. The norm takes the form:
\|Ax — b)| =
In our example the matrixes A and b are:
gn 92 93 fra
A= ucr, oo ucrg—-ucr3 5 b=] term
ucpry ucpr2 ucpr3 tepr
The looked up solution x is:
rm.
z=] rm
rm3
The problem (2) consists in minimizing norm ||Az —}|| without limitation on «. Matrix A
is created by extending matrix A by matrix E, where:
100
010
001
Matrix 5 is created by extending matrix b by three elements of large value (b4, bs, bg). In our
example the matrixes A, 6 and « are:
n 92 93 frd
ucry ucr2 ucrs term
=z ucpr, ucpra ucpr3 r tepr ea
ae 1 0 0 _ OS by J? | Tm
0 1 0 bs rms
0 0 1 be
In literature of mathematics we can find (see [6]) the statement that: minimizig the norm
of system |/Ax — 6|| is compensated by calculating matrix x by the formula:
a=(A’-A)?.A™-5. (x)
So, is the idea applied by authors in presented example. The benefits of such solution of
problem are that we can embedding optimization in simulation on System Dynamics model
by using language Dynamo. The technical details reader can find in program in Appendix A.
In each step during simulation, variable frd is changed, so the program simulats the solving
of x, by the formula (x), step by step during horizon of simulation. Such the solution x we can
named ”dynamic”. So the embedding the optimization in simulation on System Dynamics
model, has dynamical character.
Table 1.
Parameter | Final | Original | Lower Upper
value value limit limit
tchnil 15.182 20.0 0 40
tchni2 0.000. 10.0 0 40
tchni3 6.045 20.0 0 40
Initial value of objective function fob | 0.8425 - 10%
Final value of objective function fob 0.1421 - 10°
Final value of sffd 53.3649 - 107
Final value of sffrm 12.5733 » 10°
Final value of sffpr 38.3082 - 10°
Final value of lopr 65.8087 - 10°
3 The results of simulation and optimization
At first, let us present some of the results of the experiments that have been obtained from
simulation during optimization (using COSMIC and COSMOS).
Tabel 1 compiles the results of the first experiment. The number of iteration was 30.
Some of these results are presented in an illustrative form in Figures 4, 5 and 6.
Now, let us present some of the results of the experiments that we achieved from embed-
ding optimization in simulation (using Professional Dynamo 4.0). In Figures 7, 8 and 9 the
characteristics of rates: rm1, rm2, rm3 and characteristics of variables swm1, sum2, sum3
and summ are shown. It is interesting to compare the obtained value of objective function
fob (in experiments in COSMOS) and the obtained value of variable swmm (in experiments
in Dynamo). These variables are derived from two different philosophies. In one word, we
can say that: fob like as the classic objective function measures some aspects of the behavior
of the system in the whole simulation horizon. On the contrary, the variables sum1, sum2,
sum3 measure, “step by step”, the differences between the equations of balance and, surely,
the variable swmm is the summarizer of this differences in the whole horizon. To compare
the dynamics of both: fob and summ, see Figures 4 and 9.
Sh Puncrin OF FETTNe CoetaRCE SINE
EFC Srl FoUn TSH POS BFE Tous] FI toeeTe-Bara7erspha]_F2 print Eerser
Figure 4. The dynamics of the characteristics of variables sffd, sffpr, sffm and fob (expe-
riments in COSMOS)
10
model
BASIC MODEL OF BALANC!
30400 _-oD _-aay Pa oo
iebo we TaN SIMLLATEON ‘lira °° heart
eR
Figure 5. The dynamics of the characteristics of main levels of system (experiments in CO-
SMOS)
a COSMIC model
BASIC MODEL OF BALANCE
sot oot
emt we
o> 0
Sar
wot,
jo et —COO GOO ae =e,
THe EDK TNE MATION StHULATEON Fain Sveenb
Be «ua seno 4a azaz
nea MEE sas eR
es aE. “ans “Zins
Figure 6. The dynamics of the characteristics of variables rs1, rs2, and rs3 (experiments in
COSMOS)
‘TA.
see mt ese rm2 — m3
IAN A’
f [
VW, NX
AD ch cane a ee etale cenen ace
8 0 20 30 0 808970 80 90
TIME
Figure 7. The dynamics of the characteristics of variables rm1, rm2 and rm3 (experiments
in Professional Dynamo 4.0).
anog sum1(0.,10.e3) sss» sum2(0.,400.e6) === sum3(0.,20.e6)
20.¢1
7500 wo
300.¢
15.
Lo
5009 —]
200.6
10.e
2500, Fa
100.et
5000. ej
° 0 0 20 30 0 0 60 0 80 90 104
. TIME
Figure 8. The dynamics of the characteristics of variables sum1, sum2 and sum3 (experiments
in Professional Dynamo 4.0).
400.eg— Simm
300.
—
200.6
|
|
100.¢ —
Le
La
5020-3040 50607080 9800
TIME
Figure 9. The dynamics of the characteristics of variable swmm (experiments in Professional
Dynamo 4.0).
12
The economics interpretation of the results presented on Figure 7 we try to explain on
example of the value of variables rm1, rm2 and rm3 in time ”day 35”. The variables takes
the values:
rml = 8.75,
rm2 = 5.22,
rm3 = 13.74.
The first simulated equation of system R = Ax —B gives:
8.75 + 2-5.22 + 3-13.74 — 64.62 = —4.21, (J)
the second equation gives:
100 - 8.75 + 50 - 5.22 + 10 - 13.74 — 2700 = 1273.6, (1)
the third equation gives:
500 - 8.75 + 500 - 5.22 + 100 - 13.74 — 8000 = 359. (IID)
The result (J) has such the meaning that optimizing solution requires the limitation of produc-
tion (in relation to plan (forecasted) values). The result (IJ) indicates that optimized solution
gives cheaper costs of raw material (the value of ’savings” is 1273.6). The result (IZJ) has
the meaning that optimizing solution (in sense of norm ||Az — }j|), gives the little expensive
cost of production (peoples and technology, the value of 359). In summary, the optimizing
solution allows to save the value:
1273.6 — 359.0 = 914.6.
The one problem is that the limitation of production can lead to decrease of inventory and in
consequences, sometimes, can make impossible to satisfy the demand on market. But such is
the benefit of connecting the optimization with simulation on System Dynamics models, that
the simulation allows to study the influences of optimized solution on dynamics of a whole
system.
4 The concept of embedding optimization in simulation on
System Dynamics models and vice versa — Conclusions
The relations between simulation and optimization are illustrated in Figures 10 and 11.
In block diagrams in Figure 12 and 13 the most important elements of both formulations
are contained. The differencies between formulations may be easily observed.
Simulation
on System Dynamics }“—
models
td X
Optimization
Figure 10. Embedding simulation in optimization (in Coyle’s sense)
13
Optimization
a s
S; o>
ly hk
tion ic
on System Dyna
Figure 11. Embedding optimization in simulation on System Dynamics models (in the au-
thors’ understanding)
After presenting this synthetic form of the concept of embedding simulation in optimiza-
tion and vice versa, we have come up to the following conclusions:
yo
20
3°
4°
Both formulations have different possibilities and require different tools for their reali-
zation.
The actual, simulated in our experiments, version of model DYNBALANCE(3-1) has
many simplification, for example:
a) productivities of technology 1, 2, 3 can be defined like tables (not parameters) and
the optimization procedure will find the fully dynamic characteristic of optimal
solution (supplies of three raw materials in case I);
&
sources of raw materials can be defined like any dynamics characteristic (not only
*step” like in case I and II);
&
price of product (in case I) can be variable not parameter in whole horizon of
simulation;
&
weight factors in objective function (in case I) can take optional values to model
authors preferencies to the factors pf objective function (of course the preferencies
can change in horizon of simulation).
The minimizing of objective function fob in case I of model DYNBALANCE(3-1) is not
the only possibility to ask the question of optimization in problem of balance of raw
materials. The contrary we can formulate the problem in oposite way: maximizing the
profit from sale at the minimum discrepancies from plans (”"mass” plan, cost plans).
The tools for simulation experiments, chosen by authors, seems to be adequate to taking
the problem of optimization in both cases. Professional Dynamo 4.0 make possible to
optimizing in simulation on System Dynamics models, and COSMIC and COSMOS
allows to embedding simulation in optimization in Coyle’s sense.
Authors are at the begining of the way of incorporation the other methods with System
Dynamics and expects the discussion on that subject.
4
a ae
/print OF, parameters \
\and of the last iteration?
determine ordinary parameters
determine form of objective function OF
(min,max) for optimized parameters
time shift and relabelling in simulation
determine the number of iterations (ND)
¥
[determine the size of step multiplier]
v
value of current iteration: variable /TERATION
(initial value of JTERATION=0)
Vv
[ITERATION=ITERATION+1}¢
$< if ITERATION=NI >
no
v
simulate the dynamics of the system in
accordance with 2 to the end of "length"
¥
calculate the value of OF and parameters
print the value of OF and parameters at the end iteration
Figure 12. Block diagram of simulation during optimization
15
a,
\. time=0 /
eee
Vv
determine the structure of the system,
matrixes for optimized decision
rules: M, b
Vv
time shift as second block on Figure 10
t=t+dt
v
simulation of dynamics
of System Dynamics model
pee ” ss
1>H yes ( print the values “NS
7 \of chosen variables /
no ~
v
solving the matrix system: M x=b,
x=(M M) M b,x >0
Vv
simulation dynamic of the system according
to the structure of flows and time shift
v
if calculations of state
variables and rates for ¢
are finished
no
Figure 13. Block diagram of embedding optimization in simulation on System Dynamics
models
16
References
[1] COSMIC and COSMOS user manuals, The COSMIC Holding Co., Shrivenham, 1994.
[2] R. G. Coyle, System Dynamics modelling. A practical approach, Chapman & Hall, Lon-
don, 1996.
[3] R. G. Coyle, The practice of System Dynamics: milestones, lessons and ideas from
30 years experience, System Dynamics Review, 14 (1998), pp. 343-365.
[4] R. G. Coyle, Simulation by repeated optimization, Journal of the Operational Research
Society, 50 (1999), pp. 429-438.
[5] E. Kasperska, E. Mateja-Losa, D. Slota, Some Dynamics Balance of Production via
Optimization and Simulation within System Dynamics Method, in: Proc. 19th Int. Conf.
of the System Dynamics Society, ed. J. H. Hines, V. G. Diker, R. 5. Langer, J. I. Rowe,
Atlanta 2001, pp. 1-18.
[6] J. Legras, Praktyczne metody analizy matematycznej, WNT, Warszawa 1974 (translation
from French: Methodes et Technique De l’Analyse Numerique, Dunod, Paris 1971).
[7] R. Lukaszewicz, Management System Dynamics, PWN, Warsaw, 1975 (in Polish).
[8] R. Lukaszewicz, The direct form of structure models within System Dynamics, Dynamica,
2 (1976), pp. 36-43.
[9] Professional DYNAMO 4.0 for WINDOWS. Reference manual, Pugh-Roberts Associa-
tes, Cambridge, 1994.
Appendix A. Program in DYNAMO
* Balance of 3 raw materials
note
note level of raw material during transformation
note
n lmt=300
1 lmt.k=lmt.j+dt*(rrm. jk-rpr. jk)
note
note input rate to lmt (rrm)
note
rv rrm.kl=gi*prm1.k+g2*prm2.k+g3*prm3.k
ce gi=1
c g2=2
c g3=3
note
note proofed rate of ist raw material (rm1)
note proofed rate of 2nd raw material (rm2)
note proofed rate of 3rd raw material (rm3)
note
a prmi. .k,source1.k,source1.k,rm1.k)
a prm2. .k,source2.k,source2.k,rm2.k)
a prm3.k=clip(rm3.k,source3.k,source3.k,rm3.k)
a sourcel.k=step(100,0)
17
a source2.k=step(100,0)
a source3.k=step(100,0)
note
note output rate from lmt (rpr)
note
rv rpr.kl=lmt.k/tpr
note
note time of production (tpr)
note
c tpr=2
note
note level of inventory of production (lin)
note
-k=Lin. j+dt*(rpr. jk-rsl. jk)
note
note output rate from lin - rate of sale (rsl)
note
rv rsl.kl=clip(0,rd.k,0,lin.k)
note
note rate of demond (rd)
note
a rd.k=p0+p1*sin((6.28*time.k)/perd)
c p0=100
c p1i=30
c perd=52
note
note
note
a frd.k=p0+p2*sin((6.28*time.k)/perd)
c p2=40
note
note unit cost
note unit cost
note unit cost
note
¢ ucri=100
© ucr2=50
¢ ucr3=10
note
note unit cost
note unit cost
note unit cost
note
¢ ucpri=500
¢ ucpr2=500
¢ ucpr3=100
note
of
of
of
of
of
of
production
production
production
production
production
production
of 1st raw material (ucr1)
of 2nd raw material (ucr2)
of 3rd raw material (ucr3)
from ist raw material (ucpr1)
from 2nd raw material (ucpr2)
from 3rd raw material (ucpr3)
note total cost of raw material - plan (tcrm)
note total cost of production - plan (tcpr)
18
note
¢ tcrm=2700
c tcpr=8000
note
note first three row of matrix a
note
a ali.k=g1
a al2.k=g2
a al3.k=g3
a a2i.k=ucri
a a22.k=ucr2
a a23.k=ucr3
a a31.k=ucpri
a a32.k=ucpr2
a a33.k=ucpr3
note
note vector b
note
a bil.k=frd.k
a b2.k=tcrm
a b3.k=tcpr
a b4.k=8000
a b5.k=80000
a b6.k=16000
note
note vector bb=at.b
note
a bb1.k=al1.k*b1.k+a21.k*b2.k+a31.k*b3.k+b4.k
a bb2.k=a12.k*b1.k+a22.k*b2.k+a32.k*b3.k+b5.k
a bb3.k=a13.k*b1.k+a23.k*b2.k+a33.k*b3.k+b6.k
note
note matrix c=at.a
note
c1i.k=1+al1.k*al1.k+a21.k*a21.k+a31.k*a31.k
c12.k=a11.k*a12.k+a21.k*a22.k+a31.k*a32.k
c13.k=a11.k*a13.k+a21.k*a23.k+a31.k*a33.k
c21.k=c12.k
c22.k=1+a12.k*al12.k+a22.k*a22.k+a32.k*a32.k
c23.k=a12.k*a13.k+a22.k*a23.k+a32.k*a33.k
c31.k=c13.k
c32.k=c23.k
c33.k=1+a13.k*a13.k+a23.k*a23.k+a33.k*a33.k
sppppp pp p
note determinant of matrix c=at.a
note
a detc.k=-c13.k*c22.k*c31.k+c12.k*c23.k*c31.k+c13.k*c21.k*c32.k-7
c11.k*c23.k*c32.k-c12.k*c21.k*c33.k+c11.k*c22.k*c33.k
note
note matrix d=Det[c]*Inverse[c]
note
19
d11.k=c22.k*c33.k-c23.k*c32.k
d12.k=c13.k*c32.k-c12.k*c33.k
d13.k=-c13.k*c22.k+c12.k*c23.k
d21.k=-c21.k*c33.k+c31.k*c23.k
»k=-c13.k*c31.k+c11.k*c33.k
d23.k=c13.k*c21.k-c11.k*c23.k
d31.k=-c22.k*c31.k+c21.k*c32.k
d32.k=c12.k*c31.k-c11.k*c32.k
d33.k=-c12.k*c21.k+c11.k*c22.k
sop pp ppp p
a
vd
i
note rate of 1st raw material (rm1)
note rate of 2nd raw material (rm2)
note rate of 3rd raw material (rm3)
note rm=(d*bb) /Det [c]
note
a rmi.k=(bb1.k*d11.k+bb2.k*d12.k+bb3.k*d13.k)/detc.k
a rm2.k=(bb1.k*d21.k+bb2.k*d22.k+bb3.k*d23.k)/detc.k
a rm3.k=(bb1.k*d31.k+bb2.k*d32.k+bb3.k*d33.k)/detc.k
note
note
note
a bli.k=(b1.k-al1.k#rm1.k-a12.k*rm2.k-a13.k*rm3.k)**2
a b12.k=(b2.k-a21.k#rm1.k-a22.k*rm2.k-a23.k*rm3.k)**2
a b13.k=(b3.k-a31.k#rm1.k-a32.k*rm2.k-a33.k*rm3.k)**2
note a bl4.k=(b4.k-rm1.k)**2
note a b15.k=(b5.k-rm2.k)**2
note a b16.k=(b6.k-rm3.k)**2
a
a
note a functio.k=bli.k+b12.k+b13.k+b14.k+b15.k+b16.k
note a funi3.k=bl1.k+bl2.k+b13.k
note
note
1 sumi.k=sum1.j+dt*bli.j
1 sum2.k=sum2. j+dt*bl2.j
1 sum3.k=sum3. j+dt*b13.j
n
n
n
a summ.k=(sum1.k+sum2.k+sum3.k)
note
note parameters of simulation
note
spec length=104/dt=1/savper=1
save rm1,rm2,rm3,sum1,sum2,sum3,summ
20
Appendix B. Program in COSMOS
note
note
note
1mt.k=1mt. j+dt*(rrm. jk-rpr. jk)
1imt=0
rrm.kl=qi*rmi.kl+q2*rm2.k1+q3*rm3.k1
qi=1
q2=2
q3=3
rm1.kl=clip(tchn1,source1.k,source1.k,tchn1)
source1.k=step(100,0)
tchni=20
rm2.kl=clip(tchn2,source2.k,source2.k,tchn2)
source2.k=step(100,0)
tchn2=10
rm3.kl=clip(tchn3,source3.k,source3.k,tchn3)
source3.k=step(100,0)
tchn3=20
rpr.kl=lmt.k/tpr
tpr=2
lin. k=lin. j+dt*(rpr.jk-rsl. jk)
1in=300
rsl.kl=clip(0,rd.k,rd.k*dt,lin.k)
rd.k=p0+p1*sin(6.28+*time.k/perd)
p0=100
p1=30
perd=52
wit
w2=1
w3=1
scalei=1
scale2=1
scale3=1
scale4=1
sffid.k=sffd.j+dt*((rsi.jk)**INT(2))
sffd=wp1
wp1=0
sffrm.k=sffrm. j+dt*((rs2. jk) **INT(2))
sffrm=wp2
wp2=0
sffpr.k=sffpr.j+dt*((rs3.jk)**INT(2))
sffpr=wp3
wp3=0
rsi.kl=frd.k-q1*rm1.kl-q2*rm2.k1-q3*rm3.k1
rs2.kl=tcrm-ucri*rm1.kl-ucr2*rm2.kl-ucr3*rm3.k1
rs3.kl=tcpr-ucpri*rm1.kl-ucpr2*rm2.kl-ucpr3*rm3.k1
frd.k=p0+p2*SIN(6.28*time.k/perd)
p2=40
OPKHHHOBRPOBHFPOB RPA XD XDA DMM MOO PKB rH OKROPKHAOPKHROAHPHOAAAHKB HB
21
ucri=100
ucr2=50
ucr3=10
tcrm=2700
ucpri=500
ucpr2=500
ucpr3=100
tcpr=8000
cena=1000
kara=1
lopr.k=lopr. j+dt*(rlopr. jk)
lopr=0
dod. k=clip(rd.k-rsl.k1,0,rd.k-rsl.k1,0)
rlopr.kl=dod.k*cena
note =
note objective function
note ==
a fob. k=((wi*sffd.k)/scale1)+((w2*sffrm.k)/scale2)
x +((w3*sffpr.k)/scale3)+((kara*lopr.k)/scale4)
note ==
note output control sector
note =
ce dt=1
c length=104
c prtper=5
c pltper=5
print 1)rm1,rm2,rm3,fob
print 2)1lmt,lin,lopr,rlopr
print 3)sffd,sffrm,sffpr
print 4)rrm,rsl,rsi,rs2,rs3
plot fob=1,frd=2/sffd=3(0,1000)/sffpr=4(0,90E+08) /sffrm=5
note plot rmi=1,rm2=2,rm3=3
plot 1mt=1(0,10000) ,lin=2/lopr=3
note plot sffd=1,sffrm=2,sffpr=3
plot rsi=1(-60,60)/rs2=2,rs3=3
run basic model of balance
H@PBrRPAAAAAAAA AA
note =
note definition of variable
note =
Imt=(unit) level of material during transformation
rrm=(unit/week) rate of material
rpr=(unit/week) rate of material
tpr=(week) time of production
rd=(unit/week) actual rote of demand
lin=(unit) level inventory of production
sffd=(unit**2/week) sum function of fitting demand
scale1=(unit**2/week)
sffrm=($**2/week) sum function of fitting cost balance of raw material
scale2=($**2/week)
sffpr=($**2/week) sum function of fitting cost balance of production
aaaananaanaaaga
22
scale3=($**2/week)
rsi=(unit/week)
rs2=($/week)
rs3=($/week)
frd=(unit/week) rate of demand
rmi=(unit/week) source of raw material
rm2=(unit/week) source of raw material
rm3=(unit/week) source of raw material
source1=(unit/week)
source2=(unit/week)
source3=(unit/week)
tchni=(unit/week) productivity of technology 1
tchn2=(unit/week) productivity of technology 2
tchn3=(unit/week) productivity of technology 3
tcrm=($/week) total cost of raw material
tcpr=($/week) total cost of production technology and people
fob=(1) objective function
qi=(1) fraction of material 1
q2=(1) fraction of material 2
q3=(1) fraction of material 3
rsl=(unit/week)
note = = =
d perd=(week) parametr period of output
d time=(week) time within simulation
week) simulated period
anaeanananaAnRAARAAAAAAAAARAA
d p0=(unit/week) parametr step of input
d pi=(unit/week) parametr of amplitude of sinusoidal input
note d gi=(1)
note d g2=(1)
note d g3=(1)
ucr1=($/unit)
ucr2=($/unit)
ucr3=($/unit)
ucpri=($/unit)
ucpr2=($/unit)
ucpr3=($/unit)
cena=($)
kara=(1)
wi=(1)
w2=(1)
w3=(1)
wp1=(unit**2/week)
wp2=($+*2/week)
wp3=($+*2/week)
dod=(unit/week)
rlopr=($*unit/week)
lopr=($*unit)
scale4=($*unit)
Qaanananaananananan
