Введение в анализ, синтез и моделирование систем

Пример. Математическое и компьютерное моделирование подробно, поэтапно, мы рассмотрим на примере следующей простой модели производства. Итак, возьмем укрупненные этапы моделирования производства.

Этап 1. Содержательная постановка задачи

Современное производство характерно тем, что часть производимой продукции (в стоимостном выражении) возвращается в виде инвестиций (т.е. части конечной продукции, используемой для создания основных фондов производства) в производство. При этом время возврата, ввода в оборот новых фондов может быть различным для различного рода производства. Необходимо промоделировать эту ситуацию и выявить динамику изменения величины основных фондов производства (капитала).

Сложность и многообразие, слабая структурированность и плохая формализуемость основных экономических механизмов, определяющих работу предприятий, не позволяют преобразовать процедуры принятия решений в экономической системе в полностью эффективные математические модели и алгоритмы прогнозирования. Поэтому целесообразно использование простых, но гибких и надежных процедур принятия решения.

Рассмотрим одну такую простую модель социально-экономического процесса.

Этап 2. Формулировка гипотез, построение, исследование модели

Динамика изменения величины капитала определяется в нашей модели, в основном, простыми процессами производства и описывается так называемыми обобщенными коэффициентами амортизации (расхода фондов) и потока инвестиций (часть конечного продукта, используемого в единицу времени для создания основных фондов). Эти коэффициенты - относительные величины (оцениваются за единицу времени). Необходимо разработать и исследовать модель динамики основных фондов. Считаем при этом допустимость определенных гипотез, определяющих систему производства.

Пусть x(t) - величина основных фондов (капитала) в момент времени t, где 0tN. Через промежуток времени Δt она будет равна x(t+Δt). Абсолютный прирост равен Δx=x(t+Δt)-x(t). Относительный прирост будет равен x=[x(t+Δt)-x(t)]/Δt.

Примем следующие гипотезы:

1. социально-экономические условия производства достаточно хорошие и способствуют росту производства, а поток инвестиций задается в виде известной функции y(t);
2. коэффициент амортизации фондов считается неизменным и равным m, и при достаточно малом значении Δt, изменение основных фондов прямо пропорционально текущей величине капитала, т.е. dx=y(t) - mx(t).

Считая Δt0, а также учитывая определение производной, получим из предыдущего соотношения следующее математическое выражение закона изменения величины капитала - математическую модель (дифференциальное уравнение) динамики капитала:

x´(t) = y(t) - mx(t),   x(0)=х0,

где х(0) - начальное значение капитала в момент времени t=0.

Эта простейшая модель не отражает важного факта: социально-экономические ресурсы производства таковы, что между выделением инвестиций и их введением и использованием в выпуске новой продукции проходит время (лаг). Учитывая это, можно записать модель в виде

x´(t) = y(t-T)-mx(t),   x(0)=х0

Этой непрерывной, дифференциальной, динамической модели можно поставить в соответствие простую дискретную модель:

хi+1=хi +yj - mхi  ,   x0=с,   i=0, 1, 2, :, n,   0<j<n,           

где n - предельное значение момента времени при моделировании.

Дискретная модель следует из непрерывной при Δt=1, при замене производной x´(t) на относительное приращение (из определения производной, это справедливо при малых значениях Δt).

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

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

1. Ввод входных данных для моделирования: с=х(0) - начальный капитал; n - конечное время моделирования; m - коэффициент амортизации; s - единица измерения времени; y - инвестиции.
2. Вычисление xi от i=1 до i=n по рекуррентной формуле, приведенной выше.
3. Поиск стационарного состояния - такого момента времени j, 0jn, начиная с которого все х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);