Suboptimal control of the system, described by forresters stochastic modelof the world dynamics.pdf Модель мировой динамики Джея Форрестера хорошоизвестна [1]. Она была опубликована в 1971 г. и связываламежду собой в динамике 5 глобальных переменных-уровней(численность населения планеты, уровень загрязнения,фонды, или капиталовложения, природные ресурсы, долюфондов в сельском хозяйстве), 7 переменных-темпов (темпрождаемости, темп смертности, темп образования загрязне-ния, темп разложения загрязнения, темп фондообразования,темп износа фондов, темп использования природных ресур-сов) и ряд переменных-параметров (качество жизни, уро-вень питания, материальный уровень жизни, эффективностьфондов, относительную плотность населения, относитель-ную величину фондов, относительную величину фондов всельском хозяйстве, предписываемую уровнем питаниячасть фондов, долю капиталовложений в зависимости откачества жизни, относительное загрязнение и др.). Измене-ния во времени переменных-уровней описывались линей-ными разностными уравнениями первого порядка. Зависи-мости переменных-темпов и переменных-параметров другот друга и от переменных-уровней вОбъединим векторы x(k) и p(k) в единый век-тор-столбец y(k) переменных задачи:( ) ( ) , ( )y k =⎡⎣x kTp kT⎤⎦T, (2)так что( ) [ , , , , , , , , ,, , , , , , ,, , , , , ].Ty k P POL CI NR CIAF QL BR DR POLGPOLA CIG CID NRUR FR MSL ECIRCR CIRA CIR POLR CFIFR CIQR=С учетом (1) получаемy(k)=[x(k)T,ϕ(k,x(k))T]T. (3)Таким образом, вектор y(k) всех переменных мо-дели полностью определяется вектором состоянияx(k). Вектор состояния x(k) будем называть такжефазовым вектором.Введем n-вектор-функцию( , ( )) [ , ,, ,(1/ )( )]T ,f k x k BR DR POLG POLACIG CID NRURCIAFT CFIFR CIQR CIAF= − −− −⋅ −так что( ) 2 3 4 56 7 816 17 5, ( ) [ ( ) ( ), ( ) ( ),( ) ( ), ( ),(1/ )( ( ) ( ) ( ))]T ,f k x k p k p k p k p kp k p k p kp k p k x k= − −− −ƒ ⋅ −где ƒ = CIAFT - время задержки изменения частифондов в сельском хозяйстве (в годах). Как видим,вектор-функция f(k,x(k)) является функцией пере-менных-темпов, некоторых переменных-параметров ипеременных состояния. Но с учетом (1) эта функция,в конечном счете, есть сложная функция переменныхсостояния (фазовых переменных) x(k): f (k,x(k)) .Обозначив через h=DT шаг дискретизации повремени h=t(k+1)−t(k) k=1,N−1, запишем сис-тему уравнений динамики модели Форрестера в век-торной форме:x(k+1)=x(k)+h⋅f(k,x(k)) (4)с начальным условием(1) [ , , , , ]x = PI POLI CII NRI CIAFI T ,где при t=t(1)x1(1) = PI(для t(1) = 1900 г. PI=1, 65E+9 чел.,для t(1) = 1970 г. PI=3,6E+9 чел.,для t(1) = 2000 г. PI=5,3E+9 чел.),x2 (1) = POLI ,x3 (1) = CII ,x4 (1) = NRI ,x5 (1) = CIAFI .В дальнейшем при численном моделировании ианализе управляемой модели Форрестера будем ис-пользовать исходные данные из книги [1], за некото-рыми исключениями, оговоренными ниже.Примем, согласно Форрестеру [1], следующие на-чальные условия: t(1) = 1900 , t(N) = 2100 , h=DT=0,25(годы), PI=1,65E+9 (чел. в 1900 г.), POLI=0,2E+9 (ед.загрязнения), CII=0,4E+9 (ед. фондов), NRI=900E+9(ед. природных ресурсов), CIAFI=0,2 (безразмерна -начальная часть фондов в сельском хозяйстве).Нормальные значения параметров модели:BRN = 0,078 (нормальный темп рождаемости,часть/год, - примем его равным этому значению вме-сто 0.04, взятого в [1], чтобы удовлетворить значениючисленности населения PI=5,3E+9 чел. в 2000 г.),CIAFN=0,3 (нормальная часть фондов в сельском хо-зяйстве, безразмерна), CIGN=0,05 (нормальное фон-дообразование, ед. фондов/чел.год), CIDN=0,025(нормальный износ фондов, часть/год), DRN=0,028(нормальный темп смертности, часть/год), NRUN=1(нормальное потребление природных ресурсов, ед.ресурсов/год), PDN=26,5 (нормальная плотность на-селения, чел./кв.км.), POLN=1 (нормальное загрязне-ние, ед. загрязнения/чел.год), ECIRN = 1 (нормальнаяэффективность относительной величины фондов, ед.фондов/чел.), FN=1 (нормальный уровень питания,безразмерен).Приведем также значения некоторых констант:LA=135E+6 (площадь земли, кв.км.), POLS=3.6E+9(стандартное загрязнение, ед. загрязнения), ƒ == CIAFT=15 (время задержки изменения части фондов,лет), QLS=1 (стандартное качество жизни, ед. удовле-творенности), FC=1 (коэффициент питания). Осталь-ные значения констант модели и таблично заданныезависимости между переменными опускаем (они пол-ностью соответствуют работе Форрестера [1]).Уравнение состояния (4) позволяет, с учетом со-отношений (1), (2) или соотношения (3), получать извектора состояния x(k) в момент времени t(k) век-тор состояния x(k +1) в следующий момент дискрет-ного времени t(k +1) и вычислять в этот момент но-вое значение вектора переменных моделиy(k +1)=[x(k +1)T,ϕ(k +1,x(k +1))T]T. (5)Это позволяет шаг за шагом проследить изменениевсех переменных модели. Если со временем значенияпараметров модели изменяются, предусмотрено их«клиппирование» (вырезание и замена другими зна-чениями в заданные моменты времени). Тем самымобеспечивается необходимая нестационарность моде-ли, изменчивость ее параметров. Именно эта особен-ность модели отражается в первом аргументе функ-ции f (k,x(k)) , выражающем явную зависимостьуравнения состояния (4) от времени.СТОХАСТИЧЕСКАЯ МОДЕЛЬВекторное уравнение состояния мировой динами-ки (4), соответствующее модели Форрестера, не со-держит случайно изменяющихся факторов, являетсядетерминированным. Введение в модель стохастиче-ских возмущений позволило бы «проигрывать» болеереалистические траектории мировой динамики. В свя-зи с этим введем в правую часть уравнений (4) слу-чайные флуктуации переменных-темпов в виде гаус-совских «белых» (некоррелированных во времени)последовательностей случайных векторов ƒ(k) с ну-левыми математическими ожиданиями M{ƒ(k)}=0 идисперсионными матрицами D(k)=M{ƒ(k)ƒ(k)T},где M{⋅} - оператор математического ожидания:x(k+1)=x(k)+h(f(k,x(k))+ƒ(k)). (6)При стационарных стохастических возмущенияхматрицы D(k) = D остаются постоянными во време-ни. Однако в общем случае эти матрицы могут изме-няться. Структура матриц D(k) должна быть такой,чтобы с изменением шага h дискретизации по време-ни флуктуации переменных состояния модели не из-менялись. Это обеспечивается введением обратнопропорциональной зависимости дисперсий Dii(k)компонент ƒi(k) вектора ƒ(k) от шага h [2]:Dii(k) = ƒi2(k)/h, i= 1,n, (7)где ƒ1(k) - среднеквадратичное отклонение (СКО)темпа рождаемости-смертности; ƒ2(k)- СКО темпазагрязнения; ƒ3(k) - СКО темпа фондообразования;ƒ4(k)- СКО темпа потребления природных ресурсов;ƒ5(k) - СКО темпа изменения части фондов в сель-ском хозяйстве.При моделировании стохастической динамики бу-дем для простоты считать возмущения темпов ста-ционарными с постоянной диагональной дисперсион-ной матрицей2 2 2 2 2D=diag⎡⎣ƒ1/h,ƒ2/h,ƒ3/h,ƒ4/h,ƒ5/h⎤⎦при следующих ориентировочных значениях СКОтемпов:ƒ1(k) =0,5E+8; ƒ2(k)=1E+4; ƒ3(k) =1E+4;ƒ4(k)=1E+4; ƒ5(k) =1E-2.ПОСТАНОВКА ЗАДАЧИ ОПТИМАЛЬНОГОУПРАВЛЕНИЯПусть для наблюдения доступен весь вектор со-стояния x(k) в каждый данный момент дискретноговремени как в детерминированной, так и в стохасти-ческой системе. Пусть при этом наблюдения произво-дятся без ошибок. Это довольно сильное и, по-видимому, не вполне реалистичное предположение.Но для выяснения возможностей синтеза управленийв такой сложной системе, как мировая динамика, по-добное упрощение представляется оправданным.Поставим задачи синтеза оптимального управле-ния системой, описываемой детерминированной мо-делью Форрестера (4) или более общей стохастиче-ской моделью (6).Одним из основных выводов работы Форрестера[1] является тревога, вызванная постоянным ростомчисленности населения планеты, прогрессирующимзагрязнением окружающей среды, истощением при-родных ресурсов, уменьшением площадей для сель-скохозяйственной деятельности и другими явления-ми, могущими привести к возможному существенно-му снижению качества питания и уровня жизни лю-дей. Мировая динамика нестабильна, далека от равно-весия. В связи с этим возникает вопрос, нельзя ливвести в модель мировой динамики стабилизирующиефакторы, глобальные управления, способные удержи-вать ее в состоянии некоторого равновесия.Очевидно, такие управляющие факторы могутбыть введены только в правую часть уравнений дина-мики (4) или (6). В простейшем случае эти управле-ния могут линейно влиять на темпы изменения пере-менных состояния, то есть входить в правую частьуравнений состояния аддитивно:x(k+1) =x(k)+h(f(k,x(k))+ƒ(k)+Bu(k)), (8)где u(k)- вектор управляющих воздействий; B - мат-рица передачи управлений соответствующей размер-ности. Управления u(k) могут влиять на темпы ростанаселения, темпы загрязнения, темпы роста фондов,темпы использования природных ресурсов и, нако-нец, на долю капиталовложений в сельское хозяйство.Для синтеза оптимального управления необходимозадать цель управления, обычно связываемую с ми-нимизацией некоторого аддитивного функционала натраектории переменных системы и управлений:( ) ( )11 0 ,[ , ] , ( ), ( ) , ( ) minNk u xJ y u f k y k u k N y N−==ƒ +ƒ  .В качестве цели управления можно выбрать ста-билизацию на каком-либо заданном уровне одной илинескольких переменных (переменных состояния, пе-ременных-темпов или переменных-параметров). Вэтом случае можно выбрать квадратичный критерийуправления:( ) ( ) 1* *1[ , ] 1 ( 1) ( 1) ( 1)2N TkJ y u y k y R k y k y−== ƒ + − + + − +11 ,1 ( ) ( ) ( ) min2NTk u yu k Q k u k−=+ ƒ  , (9)где вектор y* определяет желаемые уровни стабили-зации соответствующих переменных. Неотрицательноопределенная матрица R(k +1) «штрафует» отклоне-ния от желаемых значений тех переменных, которыетребуется стабилизировать, обнуляя вклад в функ-ционал тех переменных, стабилизация которых нетребуется. Матрица Q(k) квадратичной формы поуправлениям должна быть положительно определен-ной (следовательно, невырожденной). Она «штрафу-ет» большие значения управлений. Чем больше диа-гональные элементы этой матрицы, тем меньше ре-сурс соответствующих компонент управления, и на-оборот.Синтез оптимального управления нелинейной сис-темой, тем более при неаналитическом задании нели-нейностей, в общем случае является крайне сложнымделом даже при таком простом критерии оптимально-сти, как квадратичный. Однако для линеаризованнойподходящим образом системы в отсутствие ограниче-ний на управление найти структуру оптимальногоуправления не составляет труда, поскольку решениезадачи оптимального управления линейной системой(как с непрерывным, так и с дискретным временем)при квадратичном критерии хорошо известно [3,4].Это задача аналитического конструирования опти-мального регулятора (АКОР). Полученное при этомоптимальное управление в k-й момент времениu(k) =u*(k) не стеснено никакими ограничениями,кроме квадратичной штрафной функции в функцио-нале (9). При наличии дополнительных ограниченийна управления типа неравенствui(k) ≤ci(k) , ci(k) ≥0, i= 1,l, (10)где l - размерность вектора управления, оптимальныеуправления соответственно ограничиваются:ui(k) =ci(k)sign(ui*(k)) , i= 1,l, (11)где sign (⋅) - знаковая функция. Заметим, что соот-ношение (11) определяет оптимальное управление вусловиях ограничений (10) только в случае скалярныхуправлений, когда l = 1. В векторном случае это со-отношение определяет лишь субоптимальное управ-ление.КВАЗИЛИНЕАРИЗАЦИЯУРАВНЕНИЙ СОСТОЯНИЯПредположим, нам удалось решить задачу (8), (9)при начальном условии x(1) , не равном x(1), но близ-ком к нему. Тогда последовательность x(1) , k= 1,N,образует оптимальную фазовую траекторию, соответ-ствующую начальному условию x(1) , а последова-тельность u(k), k=1,N−1 , - оптимальное управле-ние на этой траектории. Произведем линеаризациюсистемы (8) относительно траектории x(k), разложивфункцию f(k,x(k)) в правой части векторного урав-нения (8) в ряд Тейлора по x(k) около точки x(k):f(k,x(k))= f(k,x(k))+( , ( ))( ) ( )( ) ( ) ( ) ( )( )f k x kxk x k o xk x kx k+ − + −и ограничившись линейным приближениемx(k+1) =A(k,x(k))x(k)++h(v(k,x(k))+ƒ(k)+Bu(k)), (12)где обозначено:( ) ( , ( )), ( )( )f k x kA k x k I hx k= +; (13)( ) ( ) ( , ( )), ( ) , ( ) ( )( )f k x kv k x k f k x k x kx k= −. (14)Система (12) линейна по состояниям x(k) иуправлениям u(k). Естественно, при x(1) = x(1) этасистема в точности совпадает с нелинейной системой(8). В этом случае система (12) становится как бы ли-нейной, квазилинейной [5]. Но для такой квазилиней-ной системы (12) уже можно построить оптимальноев смысле критерия (9) управление.СТРУКТУРА АЛГОРИТМАОПТИМАЛЬНОГО УПРАВЛЕНИЯВыпишем формулы, определяющие структуру оп-тимального управления для нестационарной линейно-квадратичной задачи (9), (11) (задачи АКОР) с дис-кретным временем при точно известном начальномусловии x(1) и точных и полных наблюдениях векто-ра состояния x(k) [4, 5]:u*(k) = −Q−1(k)BTP(k+1)A(k,x(k))ƒx(k), (15)где ƒx(k)=x(k)−x* - отклонение вектора состоянияx(k)от значения x*, соответствующего желаемомузначению y* вектора y(k) переменных системы, аматрица P(k + 1) определяется решением рекуррент-ного матричного уравнения Риккати в любой из удоб-ных формP(k) = R(k)+A(k,x(k))T P(k+1)( ) ( ) I BQ 1BT P(k 1) 1 A k,x(k) + − + − ; (16')P(k)= A(k,x(k))T P(k+1)A(k,x(k))−( ) ( )1, ( ) ( 1) ( 1) ( ) A k x k T P k B BTP k B Q k−− + + + BTP(k+1)A(k,x(k))+R(k); (16'')P(k) = R(k)+A(k,x(k))T( ) ( ) P(k 1)1 BQ1BT 1 A k,x(k) + − + − − (16''')в обратном времени с граничным условием на правомконце P(N)= R(N).Если нет дополнительных ограничений на управ-ления, то u(k) =u*(k). Если имеются дополнительныеограничения вида (10), оптимальное управление кор-ректируется в соответствии с формулой (11) с учетомзамечания к этой формуле.Получив структуру алгоритма построения опти-мального управления (15), (16), перейдем к пределупри x(1)  x(1) . Тогда в выражениях (15), (16) траек-тория x(k) заменится на x(k) - оптимальную управ-ляемую траекторию. Линеаризация больше не потре-буется. Уравнения для оптимальной управляемой тра-ектории и оптимального управления примут видx(k+1) =x(k)+h(f(k,x(k))+ƒ(k)+Bu(k)); (17)u*(k) = −Q−1(k)BTP(k+1)A(k,x(k))ƒx(k); (18)ƒx(k)=x(k)−x*; (19)u(k) =u*(k), (20')или u(k) = c(k)sign(u*(k)); (20'')P(k) = R(k)+A(k,x(k))T P(k+1)( ) ( ) I BQ 1BT P(k 1) 1 A k,x(k) + − + − , (21')или P(k) = A(k,x(k))T P(k+1)A(k,x(k))−( ) ( )1, ( ) ( 1) ( 1) ( ) A k x k T P k B BTP k B Q k−− + + + BTP(k+1)A(k,x(k))+R(k), (21'')или P(k) = R(k)+A(k,x(k))T( ) ( ) P(k 1)1 BQ1BT 1 A k,x(k) + − + − − (21''')с граничными условиямиx(1) = x(1) , P(N)= R(N). (22)Получилась нелинейная двухточечная краевая за-дача (ДТКЗ). Чтобы вычислить управление (18) в те-кущий момент времени, нужно решить уравнениеРиккати (21) в обратном времени, начиная с моментаокончания процесса управления. Но матрица A(k, x(k))этого уравнения зависит от состояния системы в те-кущий момент времени, а оно может быть достигнутотолько при наличии управлений до этого момента, тоесть решений уравнения Риккати до этого момента.Следовательно, оптимальное управление, определяе-мое соотношениями (18)-(22), оказывается физическинереализуемым.ПОСТРОЕНИЕ АЛГОРИТМОВФИЗИЧЕСКИ РЕАЛИЗУЕМОГО СУБОПТИ-МАЛЬНОГО УПРАВЛЕНИЯДля реализации алгоритма оптимального управле-ния необходимо знание матрицы A(k, x(k)) в любоймомент времени до и после текущего момента, чтофизически невозможно, так как эта матрица зависитот траектории движения управляемой системы и неможет быть известной для моментов времени, буду-щих по отношению к текущему.Отсюда следует, что для построения физическиреализуемого алгоритма управления необходимо ка-ким-то образом вычислять (точно или приближенно)матрицы A(k, x(k)) на траектории x(k).Очевидно, проще всего это можно сделать, еслиуправляемая система слабо нелинейна, то есть еслиматрица A(k, x(k)) сравнительно мало изменяется нафазовой траектории x(k), в том числе на неуправляе-мой траектории. Тогда эта матрица будет мало отли-чаться от своего среднего значения на траекторииx(k), в том числе неуправляемой:A(k,x(k))≈ A, ( )11 , ( )NkA A k x kN == ƒ . (23)Рассмотрим два возможных варианта реализацииалгоритма управления на основе предположения осправедливости условия (23). Эти алгоритмы не будут,строго говоря, оптимальными, но в условиях предпо-ложения (23) они будут близкими по качеству к опти-мальному. Будем называть их субоптимальными.Алгоритм 1Предположим, что соотношение (23) действитель-но имеет место.1. Заменим в уравнении Риккати (21) матрицуA(k, x(k)) матрицей A , вычисленной по формуле (23)на неуправляемой траектории, и решим уравнение(21) в обратном времени с граничным условием (22)на правом конце. Запомним полученную последова-тельность матриц P(N),...,P(1) .2. По заданному начальному состоянию x(1) вы-числим матрицу A(1, x(1)) и, используя формулу (18),начальное управлениеu*(1) = −Q−1(1)BTP(2)A(1,x(1))ƒx(1) ,скорректировав его, если нужно, по формуле (20).3. Получив с этим управлением одношаговое ре-шение x(2) уравнения замкнутой (управляемой) сис-темы, по известному теперь состоянию x(2) вычислимматрицу A(2,x(2)) и, используя формулу (18), управ-лениеu*(2) = −Q−1(2)BTP(3)A(2,x(2))ƒx(2) ,скорректировав его, если нужно, по формуле (20).4. Продолжая этот процесс, по известному на k-мшаге состоянию x(k) вычислим матрицу A(k,x(k))и, используя формулу (18), управлениеu*(k) = −Q−1(k)BTP(k+1)A(k,x(k))ƒx(k),скорректировав его, если нужно, по формуле (20).5. Получив реализовавшееся с этим управлением всоответствии с уравнением (17) состояние x(k +1),переходим на шаг 4 с увеличением индекса k на еди-ницу.6. По достижении шага N (финального состоянияx(N) управляемого процесса) заканчиваем процессуправления.Мы получили физически реализуемый субопти-мальный алгоритм управления с обратной связью,причем матрица P(k +1), входящая в матрицу коэф-фициентов обратной связи в выражении (18), вычис-ляется заранее с использованием усредненной матри-цы A . Заменять усредненной матрицей A матрицуA(k,x(k)) в выражении (18) для управления нет не-обходимости. Она может быть вычислена в каждыйтекущий момент времени на траектории системы. Од-нако в вычислительном отношении алгоритм 1 проще,когда матрица A(k,x(k)) в выражении (18) тоже за-менена матрицей A .Таким образом, соотношения (17)-(22), опреде-ляющие структуру оптимального управления, не из-меняются для субоптимального реализуемого управ-ления, только в уравнении (21) матрица A(k,x(k))заменяется усредненной по неуправляемой траекто-рии матрицей A . Остается лишь проверить справед-ливость предположения (23).Модифицированный алгоритм 1Упрощенной модификацией алгоритма 1 являетсяалгоритм, в котором, как и в алгоритме 1, уравнениеРиккати (21) использует усредненную матрицу Aвместо матрицы A(k,x(k)) и решается в обратномвремени с условием (22) на правом конце. Но запоми-нается только матрица P0 =P(1) , и на всех шагах ал-горитма вместо матрицы P(k +1) в выражении (18)для управления используется матрица P0 . Для про-стоты матрица A(k,x(k)) в выражении (18) также за-меняется матрицей A .Алгоритм 2Другой вариант алгоритма субоптимальногоуправления на основе предположения о справедливо-сти условия (23) может быть реализован на текущейуправляемой траектории. По форме он совпадает с ал-горитмом 1, но матрицы A и P вычисляются по-иному.1. Сначала в качестве матрицы A выбирается на-чальная матрица A(1) =A(1,x(1)) и с ней в обратномвремени решается уравнение Риккати (21) с гранич-ным условием (22) на правом конце до полученияматрицы P(2) . Матрица A(1) запоминается.2. С полученной матрицей P(2) по заданному на-чальному состоянию x(1) по формулам (18), (19) на-ходится управление u(1) и, согласно уравнению (17),реализуется состояние управляемого процесса x(2) .3. По известному теперь состоянию x(2) находит-ся матрица A(2,x(2)). Усредненная по двум шагамматрицаA(2)=(A(1,x(1))+A(2,x(2)))2принимается за новую матрицу A . Матрица A(2) за-поминается.4. Снова решается уравнение Риккати (21) с мат-рицей A = A(2) на всей траектории и с граничнымусловием (22) на правом конце до получения матрицыP(3) . С полученной матрицей P(3) находится управ-ление u(2) и реализуется состояние управляемогопроцесса x(3) .5. По известному теперь состоянию x(3) находит-ся матрица A(3,x(3)) и усредненная по трем шагамматрицаA(3)=(A(1,x(1))+A(2,x(2))+A(3,x(3)))3.Матрица A(3) запоминается.6. Затем снова решается уравнение Риккати (21)с матрицей A = A(3) на всей траектории от финаль-ного момента времени до текущего с граничным ус-ловием (22) на правом конце до получения матрицыP(4) . С полученной матрицей P(4) находится љЗ†упимеющие соответственно номера 8, 9, 17, 22, равны 0,т.е. в логарифмическом масштабе равны − , поэтомуна рис. 1 они не представлены (пропущены). Эти про-пущенные значения надо понимать просто как равныенулю. На этом же рисунке приведены матрицыA ∓ S в том же представлении.Как видно из сравнения табл. 1 и 2 и представлен-ных на рис. 1 массивов, элементы матрицы A(k, x(k))относительно незначительно флуктуируют на траек-тории движения (эволюции) системы. Это дает осно-вание считать систему слабо нелинейной и предполо-жение (23) справедливым.1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 2510-2010-1510-1010-51001051010Номер элемента матрицы A|среднее A| - СКО A|среднее A||среднее A| + СКО AРис. 1. Изменчивость матрицы A на траекторииРассмотрим в качестве примера реализацию и де-терминированной и стохастической модели мировойдинамики Форрестера в отсутствие управления и вслучае использования модифицированного алгоритмасубоптимального ограниченного управления, в кото-ром управления производятся только путем перерас-пределения капиталовложений (фондов) в сельскоехозяйство, а критерием управления является стабили-зация численности населения планеты на желаемоммировым сообществом уровне.Возьмем в качестве желаемого значения числен-ности населения x1* =5E+9 (5 млрд) человек при на-чальном значении в 1900 г. x1(1)=1, 65E+9 человек.Желаемых значений остальных переменных устанав-ливать не будем. Тогда вектор рассогласования ƒx(k),определяемый формулой (19) и входящий в выраже-ние (18) для оптимального управления, примет вид*ƒx(k)=[x1(k)−x1,0, 0, 0,0]T .Матрица R(k) в этом случае должна иметь толькоодин отличный от нуля элемент - первый диагональ-ный. Примем R(k) = diag[1,0,0,0,0] k=1,N.В качестве управления u возьмем скалярную пе-ременную - темп изменения доли капиталовложенийв сельское хозяйство. В этом случае матрица B при-мет видB = [0,0,0,0,1]T .Примем Q(k)=3E+19 k=1,N−1.Ограничим управление в соответствии с формулой(20''), выбрав c(k) =0,03 k=1,N−1.Стохастические возмущения в рассматриваемомпримере будем вводить только в виде случайных не-зависимых во времени флуктуаций темпа рождаемо-сти-смертности, имеющих нормальное распределениес нулевым средним и стандартным отклонением(СКО) ƒ1(k) =0,5E+8. Тогда дисперсионная матрицастохастических возмущений в системе примет вид2D=diag ⎡⎣ƒ1 /h, 0, 0, 0, 0⎤⎦ .Нетрудно проверить, что линеаризованная система(12) при выбранной матрице B = [0,0,0,0,1]T , а следо-вательно, и нелинейная система (17), вполне управ-ляемы [3], поскольку матрица управляемости системы(12) с усредненной матрицей динамики A имеет пол-ный ранг n = 5 и матрица W(1,N) не вырождена:rank[B,AB,A2B,A3B,A4B]=5, detW(1,N)  0 ,11(1, ) ( , 1) ( , 1)NT TkW N N k BB N k−== ƒƒ + ƒ + ,ƒ(N,k+1) = A(N−1,x(N−1))⋅⋅⋅A(k +1,x(k +1)).Результаты численного моделирования на двух-сотлетнем интервале времени (с 1900 по 2100 г.) не-управляемой и управляемой систем мировой динами-ки Форрестера в отсутствие и при наличии стохасти-ческих возмущений приведены на рис. 2 - 12. Естест-венно, в силу ограниченности места мы приводим ходтолько некоторых переменных модели.1900 1920 1940 1960 1980 2000 2020 2040 2060 2080 2100-0.00500.0050.010.0150.020.0250.030.035Годыu - Субоптимальное ограниченное управлениедетерминированная системастохастическая системаРис. 2. Субоптимальное ограниченное управлениеНа рис. 2 приведен ход субоптимальных ограни-ченных управлений с обратной связью для детерми-нированной и стохастической системы. Управлениепостроено по модифицированному алгоритму 1.Видно, что на начальном (почти 20-летнем) интерва-ле времени управление принимает максимально воз-можное значение, т.е. перераспределение капиталовло-жений в сельское хозяйство производится с максималь-но возможной скоростью. Затем темп этого перераспре-деления резко снижается, даже на некоторое время ста-новится отрицательным, а затем стабилизируется насравнительно невысоком положительном уровне.Такое поведение управления легко объясняетсятребованием критерия управления стабилизироватьчисленность населения на уровне 5 млрд человек.Действительно, как видно из рис. 3, в начальныймомент времени (1900 г.) численность населения(1,65 млрд человек) еще далека от желаемой (5 млрдчеловек). Для увеличения численности населения помодели Форрестера требуется существенное увеличе-ние производства продуктов питания, что и потребо-вало в начальном периоде быстрейшего увеличениядоли капиталовложений (фондов) в сельское хозяйст-во. К концу этого периода население почти достигает,как видно из рис.3, желаемого значения численности,и доля капиталовложений в сельское хозяйство долж-на быть резко снижена, чтобы не допустить дальней-шего роста населения.1900 1920 1940 1960 1980 2000 2020 2040 2060 2080 21001234567 x 109ГодыP - Населениедетерминированная неуправляемаясреднее - СКОсреднее + СКОстохастическая неуправляемаядетерминированная управляемаястохастическая управляемаяРис. 3. Динамика численности населенияНа приведенных здесь же кривых роста численно-сти населения в неуправляемой системе стабилизациянаселения в рассматриваемом интервале времени(1900 - 2100 гг.) не наступает (возможно, стабилиза-ция произойдет значительно позже из-за сопровож-дающего рост населения истощения ресурсов).На рис. 4 представлена относительная (в расчетена душу населения) величина фондовления). По достижении некоторого минимального длясистемы значения уровень питания стабилизируется.1900 1920 1940 1960 1980 2000 2020 2040 2060 2080 21000.60.70.80.911.11.21.31.4ГодыFR - Уровень питаниядетерминированная неуправляемаястохастическая неуправляемаядетерминированная управляемаястохастическая управляемаяРис. 7. Относительный уровень питанияОчевидно, перераспределение капиталовложенийв сельское хозяйство должно привести к уменьшениюдоли фондов в другие сферы человеческой деятельно-сти, в том числе в материальную сферу. Поэтому ма-териальный уровень жизни населения должен сни-зиться, тем более что капиталовложения в материаль-ную сферу резко снижаются в начальный период ин-тервала управления, затрудняя поддержание и разви-тие материальной сферы.На рис.8 приведены кривые хода материальногоуровня жизни при наличии и отсутствии управления.Материальный уровень жизни в модели Форрестераесть безразмерная величина, которая описывает сте-пень изменения эффективности относительной вели-чины фондов на душу населения в сравнении с неко-торым стандартом - фиксированным значением.1900 1920 1940 1960 1980 2000 2020 2040 2060 2080 2100-0.4-0.3-0.2-0.100.10.20.3ГодыMSL - Материальный уровеньНа рис.11 представлены кривые изменения качест-ва жизни при наличии и в отсутствие управления.1900 1920 1940 1960 1980 2000 2020 2040 2060 2080 21000.150.20.250.30.350.40.450.50.550.60.65ГодыQL - Качество жизнидетерминированная неуправляемаястохастическая неуправляемаядетерминированная управляемаястохастическая управляемаяРис. 11. Качество жизниКачество жизни в модели Форрестера использует-ся как мера функционирования мировой системы. Ка-чество жизни определяется как степень удовлетво-ренности жизнью. Она зависит от многих факторов:уровня питания (значительно больше, чем от матери-ального уровня, когда уровень питания низок), мате-риального уровня (когда уровень питания высок, аматериальный уровень низок), плотности населения(если она слишком велика), загрязнения окружающейсреды (если оно ощутимо влияет на условия жизни).Эти и некоторые другие факторы достаточно полноучитываются в модели Форрестера [1].На рис.11 видно, что в неуправляемой системе ка-чество жизни постепенно снижается, стабилизируясь,в конце концов, на некотором не очень высокомуровне. И объяснить это можно увеличением числен-ности и соответственно плотности населения, увели-чением уровня загрязнения среды, истощением ресур-сов и т.д. В управляемой же системе, вопреки ожида-ниям, связанным со стабилизацией численности насе-ления, дело с качеством жизни обстоит еще хуже: ка-чество жизни резко падает в период стимуляции ростанаселения и устанавливается со стабилизацией чис-ленности населения на еще более низком уровне, чемв неуправляемой системе.Объяснить в какой-то мере это парадоксальное яв-ление можно с помощью рис.12, показывающего ходво времени изменения доли капиталовложений в за-висимости от качества жизни. Как видно из сравнениярис. 3, 11 и 12, при стабилизации численности насе-ления в управляемой системе путем резкого перерас-пределения фондов в сельское хозяйство уровень ка-питаловложений в остальные сферы жизни резко сни-зился, что привело к быстрому снижению качестважизни. А это, в свою очередь, привело к снижениюдоли капиталовложений, направляемых наповышение качества жизни. И хотя численность насе-ления в управляемой системе сравнительно быстропришла в желаемое стационарное состояние, качествожизни стабилизировалось на сравнительно низкомуровне, а это стабилизировало также на сравнительнонизком уровне долю капиталовложений, направляе-мых на повышение качества жизни, что опять-такистабилизирует на сравнительно низком уровне каче-ство жизни.1900 1920 1940 1960 1980 2000 2020 2040 2060 2080 21000.730.740.750.760.770.780.790.80.81ГодыCIQR - Доля капиталовложений в зав. от кач. жизнидетерминированная неуправляемаястохастическая неуправляемаядетерминированная управляемаястохастическая управляемаяРис. 12. Доля капиталовложений в зависимостиот качества жизниОднако, если бы каким-нибудь образом, например,путем технического прогресса, удалось повысить эф-фективность капиталовложений, направляемых на по-вышение качества жизни, парадокс снижения качест-ва жизни при стабилизации численности населениямог бы легко разрешиться. Следовательно, парадоксобусловлен несовершенством модели, тем обстоя-тельством, что в рассматриваемом примере все пара-метры-константы модели действительно остаются по-стоянными, тогда как в реальности они могут сущест-венно изменяться (и изменяются!) благодаря техниче-скому прогрессу. В модели просто отсутствуют пере-менные, отвечающие за технический прогресс, хотя вмодели Форрестера предусмотрена возможностьскачкообразного изменения параметров модели в лю-бой заданный момент времени с помощью «клиппи-рующих» функций.В заключение заметим, что использование другихформ алгоритма субоптимального управления (алго-ритма 1 или 2) приводит практически к тем же резуль-татам и существенно не меняет выводов, полученных врезультате численного исследования управляемой де-терминированной или стохастической модели мировойдинамики Форрестера. И главный вывод состоит в том,что модель Форрестера управляема, то есть допускаетстабилизацию вектора состояния и, следовательно,всех других переменных модели при наличии доста-точных ресурсов управления.

Скачать электронную версию публикации