Этап 3. Построение алгоритма и программы моделирования

Возьмем для простоты режим моделирования, когда m, c - известны и постоянны, y - увеличивается на каждый следующий момент времени на 1%, а также рассмотрим наиболее простой алгоритм моделирования в укрупненных шагах.

1. Ввод входных данных для моделирования: с=х(0) - начальный капитал; n - конечное время моделирования; m - коэффициент амортизации; s - единица измерения времени; y - инвестиции.
2. Вычисление x_i от i=1 до i=n по рекуррентной формуле, приведенной выше.
3. Поиск стационарного состояния - такого момента времени j, 0 j n, начиная с которого все х_j, х_{j +1}, :, х_n постоянны или изменяются на малую допустимую величину ε >0.
4. Выдача результатов моделирования и, по желанию пользователя, графика.

Алгоритм, записанный на учебном алгоритмическом языке, имеет вид

алг Производство (арг вещ m, c, n, рез вещ таб х[1:366], лит p, q);дано | производство с основными фондами, изменяющимися по закону: | х[i+1]=х[i]+y-mx[i], x[0]=c, i=0, 1, 2, :, n, 0<j<n, | t=i*h, h=1 - шаг по времени (день), | i - текущий момент времени, | m - коэффициент амортизации, | х[0]=с - заданная начальная величина капитала, | y - увеличиваемая на 1% каждый раз величина инвестицийнадо | промоделировать динамику основных фондов, т.е. выяснить: | 1) чему они равны на момент времени n; | 2) наступает ли гибель предприятия, т.е. обращается ли капитал | (основные фонды) в нуль при некотором t, и найти это t; | 3) наступает ли ситуация, когда капитал стабилизируетсянач | начало тела алгоритма | описание типов переменных цел i, | i - переменная цикла прогноза (текущее время) j, | j - задаваемая величина лага k, | k - момент гибели предприятия (если есть) y | y - величина инвестиций, увеличиваемая по закону y:=1.01*y ввод (m, n, c, y) | ввод исходных данных p:='предприятие не гибнет' | задаем начальное значение s q:='капитал не стационарен' | задаем начальное значение q х[0]:=с | начальное значение капитала (не нулевое) i:=0 | задаем начальный момент времени моделирования нц пoка (i<=n) и (х[i]>0) | заголовок цикла прогноза капитала | тело цикла прогноза капитала х[i+1]=х[i]+y-mx[i] | вычисление прибыли в следующий момент y:=1.01*y | и увеличиваем на 1% - для следующего момента если х[i+1]<=0 | проверка гибели то | если гибнет, - выполняется блок погибшего предприятия p:="предприятие гибнет" | заменяем значение s k:=i-1 | и фиксируем время гибели нц для j от k до n | цикла вычисления всех x[j]=0 | остальных, нулевых значений прибыли кц | конец блок обработки погибшего предприятия если х[i+1]=х[i] | проверка стационарности прибыли то q:="капитал стационарен" | заменяем старое значение q кцкон.

Приведем программу на Паскале для имитационного моделирования (программа реализована для функции типа y=at+b, где a, b - коэффициенты потока инвестиций; структурированность и интерфейс программы "принесены в жертву" компактности, простоте и понятности программы).

PROGRAM MODFOND; {Исходные данные находятся в файле in.dat текущего каталога}{Результаты записываются в файл out.dat текущего каталога}Uses Crt, Graph, Textwin;Type Vector = Array[0..2000] of Real; Mas = Array[0..2000] of LongInt;Var Time, Lag, t, dv, mv, i, yi, p :Integer; tmax, tmin :LongInt; a, b, m, X0, maxx, minx, aa, bb, cc, sx, tk :Real; x :Vector; ax, ay :Mas; ch :Char; f1, f2 :Text;{-------------------------------------------------------------------------------------------}Procedure InputKeyboard; { Ввод с клавиатуры }Begin OpenWindow(10,5,70,20,' Ввод данных ',14,4); ClrScr; WriteLn; WriteLn('Введите время Т прогнозирования системы:'); Repeat Writeln('Для удобства построения графика введите Т не меньше 2'); Write('Т='); ReadLn(Time); until Time>=2; WriteLn('Введите лаг:'); Repeat Write('Лаг должен быть строго меньше Т - '); ReadLn(Lag); until Lag<Time; WriteLn('Введите коэффициенты для вычисления потока инвестиций'); Write('Введите a>0: a= '); ReadLn(a); Write('Введите b>0: b= '); ReadLn(b); Repeat Write('Введите коэффициент амортизации ( 0<M<1 ) - '); Readln(m); until (m<1) and (m>0); Write('Введите значение фондов в начальный момент - '); Readln(X0); CloseWindow;end;{-------------------------------------------------------------------------------------------}Procedure InputFile; { Ввод из файла }Begin Assign(f1,'in.dat'); Reset(f1); Readln(f1,time,lag,a,b,m,X0); Close(f1);End;{-------------------------------------------------------------------------------------------}Procedure OutputFile; { Запись результата работы в файл }Begin Assign(f2,'out.dat'); Rewrite(f2); WriteLn(f2,' Результаты моделирования:'); WriteLn(f2,'Значение фондов в заданное время Т = ',x[time]:4:2); WriteLn(f2,'Максимальное значение фондов = ',maxx:4:2); Write(f2,'Минимальное значение фондов = ',minx:4:2); Close(f2);End;{------------------------------------------------------------------------------------------}Procedure InputRnd; { Ввод случайными числами }Begin Randomize; Repeat Time:=Random(90); until Time>=2; Repeat Lag:=Random(80); until Lag<Time; a:=Random(10); b:=Random(10); m:=Random; X0:=Random(200);End;{------------------------------------------------------------------------------------------}Procedure OutputScreen; { Вывод данных на экран }Begin OpenWindow(10,5,70,20,' Вывод данных: ',4,3); WriteLn; WriteLn(' Данный набор входных параметров:'); WriteLn; WriteLn(' Время Т - ',time); WriteLn(' Лаг - ',lag); WriteLn; WriteLn('Коэффициенты потока инвестиций:'); WriteLn; WriteLn(' a - ',a:4:2); WriteLn(' b - ',b:4:2); WriteLn; WriteLn('Эмпирический коэффициент амортизации - ',m:4:3); Write('Состояние фондов в начальный момент - ',X0:4:2); ReadKey; CloseWindow;End;{-------------------------------------------------------------------------------------------}Procedure Worker; { Рабочая процедура }Var yt :real;Begin x[0]:=X0; For t:=1 to Time do Begin If t<Lag+1 then yt:=0 else yt:=a*(t-1-Lag)+b; x[t]:=yt+(1-m)*x[t-1]; End; maxx:=x[0]; minx:=x[0]; tmin:=0; tmax:=0; For t:=1 to Time do If x[t]>maxx then begin maxx:=x[t]; tmax:=t; end else if x[t]<minx then begin minx:=x[t]; tmin:=t; end; OpenWindow(10,5,70,13,' Результат работы модели: ',14,7); ClrScr; WriteLn; WriteLn('Значение фондов в заданное время Т = ',x[time]:4:2); If tmin<>0 then WriteLn(' Величина фондов возрастает с ',tmin,' до ',tmax); WriteLn(' Максимальное значение фондов = ',maxx:4:2); Write(' Минимальное значение фондов = ',minx:4:2); ReadKey; CloseWindow;End;{---------------------------------------------------------------------------------------------}Procedure Mas_OX; { Масштабирование по оси ОХ }Var st :String;Begin p:=1; While Time>p*24 do inc(p); For i:=1 to 24 do Begin Str(p*i,st); OutTextXY(65+20*i,420,st) End; For t:=0 to Time do ax[t]:=70+round(20*t/p);End;{-------------------------------------------------------------------------------------------}Procedure Mas_OY; { Масштабирование по оси ОУ }Var st :String; k, r :Integer;Begin If maxx>16 then Begin k:=1; While maxx>k*16 do inc(k); For i:=1 to 16 do Begin Str(k*i,st);OutTextXY(35,407-20*i,st);End; tk:=k; End else Begin r:=1; While (maxx<=16/r) and (r<16) do inc®; dec®; For i:=1 to (trunc(16/r-0.1)+1) do Begin Str(i,st); OutTextXY(35,407-0*r*i,st) End; tk:=1/r; End; For t:=0 to Time do ay[t]:=410-round(20*x[t]/tk);End;{----------------------------------------------------------------------------------------------}Procedure Ipol(x1,y1,x2,y2,x3,y3:Real); {Процедура интерполяции}Var d1, da, db, dc :Real;Begin d1:=x1*x1*(x2-x3)+x2*x2*(x3-x1)+x3*x3*(x1-x2); da:=y1*(x2-x3)+y2*(x3-x1)+y3*(x1-x2); db:=x1*x1*(y2-y3)+x2*x2*(y3-y1)+x3*x3*(y1-y2); dc:=x1*x1*(x2*y3-y2*x3)+x2*x2*(x3*y1-y3*x1)+x3*x3*(x1*y2-y1*x2); aa:=da/d1; bb:=db/d1; cc:=dc/d1;End;{--------------------------------------------------------------------------------------------}Procedure Graf; { Построение графика }Begin dv:=detect; InitGraph(dv,mv,''); SetBkColor(7); SetColor(6); Rectangle(30,40,600,450); Line(600,60,620,60); Line(620,60,620,470); Line(50,450,50,470); Line(50,470,620,470); SetFillStyle(1,1); FloodFill(610,450,6); SetFillStyle(1,15); FloodFill(100,100,6); SetColor(5); Circle(70,410,2); Line(70,410,70,50); Line(70,410,590,410); { оси ОХ и ОУ } OutTextXY(587,407,'>'); OutTextXY(67,47,'^'); OutTextXY(57,415,'0'); OutTextXY(80,45,'X(T) - (Величина основных фондов производства)'); OutTextXY(590,415,'T'); OutTextXY(540,430,'(Время)'); SetColor(2); For i:=1 to 16 do Line(67,70+20*i,70,70+20*i); For i:=1 to 24 do Line(70+20*i,410,70+20*i,413); Mas_OX; Mas_OY; For t:=0 to time do Вegin SetColor(Blue); Circle(ax[t],ay[t],2); SetFillStyle(SolidFill,Red); FloodFill(ax[t],ay[t],Blue); End; SetColor(Red); SetLineStyle(3,1,1); Line(70,ay[time],ax[time],ay[time]); Line(ax[time],ay[time],ax[time],410); Ipol(0,x[0],1,x[1],2,x[2]); For i:=ax[0] to ax[2] do Begin sx:=p*(i-70)/20; yi:=410-round(20*(aa*sx*sx+bb*sx+cc)/tk); SetColor(Red); Circle(i,yi,1); End; For t:=1 to Time-2 do Begin Ipol(t,x[t],t+1,x[t+1],t+2,x[t+2]); For i:=ax[t+1] to ax[t+2]do Begin sx:=p*(i-70)/20; yi:=410-round(20*(aa*sx*sx+bb*sx+cc)/tk); SetColor(Red); Circle(i,yi,1); End; End; ReadKey; CloseGraph;End;{-------------------------------------------------------------------------------------------}BeginWhile true do Begin ClrScr; TextBackGround(2); Window(1,1,80,25); ClrScr; OpenWindow(30,22,50,24,' Нажмите клавишу: ',4,1); OpenWindow(5,5,75,16,' Динамика фондов производства ',14,5); ClrScr; WriteLn; WriteLn(' Пусть х(t) - основные фонды в момент времени t, y(t) -'); WriteLn(' инвестиции, m - коэффициент амортизации фондов.'); WriteLn(' Модель динамики основных фондов (L - лаг):'); Write(' x`(t) = y(t-L) - mx(t), где х(0) = Хо, y(t)=at+b, ( a,b>0 ).'); ReadKey; CloseWindow; OpenWindow(15,10,65,17,' Выбирите вариант входа-выхода: ',15,0); ClrScr; WriteLn; WriteLn(' С клавиатуры - <1>'); WriteLn(' Из файла - <2>'); WriteLn(' Случайными числами - <3>'); WriteLn(' Выход - <Esc>'); ch:=ReadKey; Сase ch of #49: InputKeyboard; #50: Вegin InputFile; OutputScreen; Еnd; #51: Вegin InputRnd; OutputScreen; End; #27: Halt(1); End; CloseWindow; Worker; OutputFile; OpenWindow(22,10,58,14,'',15,5); ClrScr; WriteLn; Write('Для просмотра графика нажмите ввод'); ch:=ReadKey; If ch=#13 then begin Graf; RestoreCrtMode; end; CloseWindow; TextBackGround(15); Window(1,1,80,25); ClrScr; OpenWindow(15,10,65,16,'',15,6); ClrScr; WriteLn; WriteLn(' Хотите еще моделировать ?'); WriteLn; WriteLn('Для выхода нажмите - < Esc >'); WriteLn('Для продолжения нажмите любую другую клавишу'); ch:=ReadKey; If ch=#27 then Halt(1); CloseWindow; End; ClrScr; TextBackGround(0);End.