Adaptive Inventory Control with Restrictions andTransport Delays .pdf Теория следящих систем применяется в задачах управления запасами в работах [1 - 4]. В настоящей работе для определения поставок на склад используются метод локально-оптимального слежения [5] и метод адаптации и прогноза на основе калмановской фильтрации и экстраполяции. Управление запасами осуществляется с учетом ограничений на транспортные средства и с учетом запаздываний, при этом осуществляется минимизации издержек на хранение товара.1. Управление запасами при точном измерении спросаПусть модель склада описывается дискретным уравнением:x(k +1) = Ax(k) + Bu(k - h) - s(k), x(0) = x0, u{j) = y(j), j = -h,-h +1,...,-1,(1)где x(k) е Rn - вектор количества продукта на складе в k-й такт (x,(k) - количество товаров i-й номенклатуры), u(k - h) е Rm - вектор заказа на поставки (и, - количество заказанного товара i-й номенклатуры), h - количество тактов транспортного запаздывания, s(k) е Rn - вектор спроса в k-м такте (s,(k) - спрос на товар i-й номенклатуры), x0 и \|/(/) (/'= -h, -h+1,..., -1) - заданные векторы. Матрицы А и В определяются характеристиками и структурой склада. Оптимизируемый локальный критерий имеет видI(k) = (x(k +1) - z)1C(x(k +1) - z) + u1 (k - h)Du{k - h),(2)где C > 0, D > 0 - весовые матрицы, z - заданный отслеживаемый вектор, «Т» -операция транспонирования. В этом разделе будем предполагать, что все компоненты вектора x(k) и s(k) измеряются точно. Вычислим значения критерия (2):I (k) = u т (k - h)( BTCB + D)u(k - h) + u T (k - h)B TC (Ax(k) - s(k) - z) +Оптимальное управление определим из условияdl (k) _ndu(k - h)Тогда, в силу (3), получим уравнение(3)(B CB + D)u{k - h) + B C(Ax(k) - s(k) - z) = 0 .(4)Выражая u(k - h) из (4), управление будет следующим:u(k - h) = -(BTCB + D)-1 BTC (Ax(k) - s(k) - z). (5) Учитывая (1), имеем следующие равенстваx(k) = Ax(k -1) + Bm (k - h -1) - s (k -1), x(k -1) = Ax(k - 2) + Bm(k - h - 2) - s(k - 2),x(k - h +1) = Ax(k - h) + Bm(k - 2h) - s(k - h) .(6) Тогда, учитывая (6), локально-оптимальное управление (5) представляется в виде u(k - h) = -(BTCB + D)-1 BTC (Ah+1 x(k - h) +AABu(k - h - i) - £ AAs(k - i) - z). (7)i=1 i=0Отметим, что управление (7) формируется в момент времени k - h и для его реализации необходимо знать состояние склада x(k - h), спроса s(k - h) и прошлые значения управлений u(k - h - i), а также необходимо осуществлять прогноз спроса для моментов времени k, k - 1,^,k - h + 1.2. Управление запасами при косвенных наблюденияхВ этом случае дополнительно введем модель спроса:s(k +1) = Rs(k) + r + q(k), s(0) = s0 ,(8)где R - постоянная матрица размерности n x n, r - постоянный вектор, q(k) - случайное возмущение. Предполагается, что имеются косвенные наблюдения за вектором спроса:w(k) = Hs(k) + т(к) ,(9) где w(k) ё RMl - вектор наблюдений; H- wixn-матрица; т (k) - случайные ошибкинаблюдений; q(k), x(k) - независимые гауссовские случайные последовательности с характеристикамиM{q(k)} = 0, M{T(k)} = 0, M{q(k) qT(/')} = ßSiy, M(x(k) тТ( j)}= 75„,(10) где M{} - математическое ожидание.Так как в постановке задачи этого раздела вектор s(k) контролируется с ошибками (см.(9)), то управление (7) можно, используя оцениватели, реализовать в видеu (k - h) = -(BT CB + Dy1 BT C(Ah+1x(k - h) + +^ AA Bu* (k - h - i) - AA sf (k - h) - £ AA sp (k - i) - z), (11)i=1 i=0где sf (k - h) - оценка фильтрации, которая определяется с помощью алгоритма оптимальной калмановской фильтрации [6]:sf (k - h) = Rsf (k - h -1) + r + Kf (k - h)[w(k - h) - H (Rsf (k - h -1) + r)],Sf (0) = So ;(12)Kf (k - h) = P(k - h / k - h -1) HT (HP(k - h / k - h -1) HT + T )-1 ;(13)P(k-h/k-h -1) = RP(k-h-1)RT + Q ;(14)P(k - h) = (E - Kf (k - h)H)P{k - h / k - h -1), P(0) = P0 ,(15)где Е - единичная матрица соответствующей размерности. Фильтр (12) использует информацию, поступившую из канала измерений в момент k - h. В (11) требуется вычислять также оценки и в моменты большие, чем k - h (оценки прогноза), поэтому здесь необходимо воспользоваться экстраполятором, который позволит вычислить оценку спроса с прогнозом на 1 такт sp (k - h +1):sp (k - h +1) = Rsp (k - h) + r + Kp (k - h)(w(k - h) - Hsp (k - h)),Sp (0) = So ;(16) Kp (k - h) = RPp (k - h) HT (HPp (k - h) HT + T )-1 ;(17)Pp (k-h+1)=(R-Kp (k-h)H)Pp (k-h)(R-Kp (k-h)H)T +Q+Kp (k-h)TKpT(k-h) ,Pp (0) = Po ,(18) а оценки прогнозов для тактов j = 2,..., h-1 определятся по формуламSp (k - h + j) = RSp (k - h + j -1) + r .(19)3. Адаптация в задаче управления запасамиВ этом случае будем предполагать, что модель спроса содержит неизвестные параметры:s(k +1) = R(Q)s(k) + r(9) + q(k) ,(20) где 9 - неизвестный постоянный вектор, Л(9) - матрица, линейно зависящая от компонент вектора 9 , r(9) - вектор, линейно зависящий от 9 . Предполагается, что имеется следующая априорная информация о векторе 9 :M{9} = 9о, M{(9-9с)(9-9с )т} = P6o.Адаптация закона управления осуществляется по принципу разделения, который в нашем случае состоит из следующих этапов:•·построение закона управления по локальному критерию в предположении, что параметры 9 в модели (20) известны точно (закон управления имеет вид (11));•·оценивание неизвестных параметров (идентификация);•·формирование закона адаптивного управления ua(k) вида (11), которое осуществляется путем замены в (12), (14), (16) - (19) точных значений параметров на их оценки.Так как параметры являются неизвестными константами, то их динамическая модель имеет вид8(*) = 8(* -1) .(21)Для идентификации будем использовать оптимальный дискретный фильтр Калмана. Учитывая модель косвенных наблюдений (9), получим модель измерений в текущий момент времени k - h в следующем виде:w(k - И) = Hs(k - И) + т(к - h) .(2)Тогда, учитывая (8), имеемw(k - И) = HR(9)s(k - h -1) + Hr(9) + Hq(k - h -1) + т(к - h) .(23)Представим (23) в линейном относительно 9 векторно-матричном виде:w(k - И) = G(s(k - h -1))9 + g(s(k - h -1)) + Hq(k - h -1) + x(k - h) ,(24)где G - матрица и g - вектор.Учитывая (21), (24), рекуррентные уравнения, осуществляющие оценивание неизвестных параметров 9 , представим в виде9(k - h) = 9(k - h -1) + Ke (к - h)(w(k - h) - G9(k - h -1) - g), 9(0) = 90 ;(25)Ke (k - h) = Pe (k - h)GT (GP (k - h)(3T + HQHT + T )-1 ;(26)Pe (k - h +1) = (E - Ke (k - h)G)Pe (k - h), Pe (0) = P6o .(27)где G = G(sf (k - h -1)), g = g(sf (k - h -1)).Уравнения, используемые для вычисления оценок и прогноза (12) - (19), в этом случае будут иметь следующий вид:sf (k - h) = Rsf (k - h -1) + r + Kf (k - h)[w(k - h) - H (Rsf (k - h -1) + r)],Sf (0) = So ;(28)Kf (k - h) = P(k - h / k - h -1) HT (HP(k - h / k - h -1) HT + T )-1 ;(29)P(k-h/k-h -1) = RP(k-h-1)RT + Q ;(30)P(k - h) = (E - Kf (k - h)H)P(k - h / k - h -1), P(0) = P0 .(31)sp (k - h +1) = Ä?p (k - h) + r + Kp (k -- h) - Hsp (k - h)), sp (0) = so ;(32)Kp (k - h) = RPp (k - h) HT (HPp (k - h) HT + T )-1 ;(33 Pp (k -h+1) = (RR - Kp (k - h)H)Ppr (k - h)(RR -Kp (k - h)H)T + ß+Kp (k - h)TKpT (k - h),Pp (0) = Po .(34) sp(k -h + j) = RR?p(k -h + j-1) + r , j = 2,..., (35) где RR = Ä(9(k - h)), r = r(9(k - h)).4. Учет транспортных ограниченийВектор поставок будем определять с учетом транспортных ограничений следующего вида (для поставок используется одно транспортное средство грузоподъемности Gmax):иа (k), если Gmin < G(wa (k)) < Gmax,и, (k) = | 0,если G(иа (k)) < G^, (36),«a (k)/ a(k)' еСЛИG(ua (k)) > Gmax nЗдесь G(ua (k)) = X /Лма,г (k) - вес перевозимого груза (pt - вес единицы товара i-й номенклатуры); ua (k) - вектор поставок, построенный на основе принципа разделения; Gmm = KrGmax (Kr - коэффициент использования грузоподъемности транспортного средства); a(k) = Gmax / G(ua (k)).5. Минимизация издержек на хранение товаровИздержки на хранение товара на заданном скользящем временном отрезке [k, k+T] определяются по формулеnk+TЛзд (z, k) = ZZ xi (t)Ci .(37)i=1 t=kПри этом должны выполняться ограниченияX(k) > wt, Vk e [k, k + T], i = 1 n,(38)где Cj - стоимость хранения единицы товара i-й номенклатуры в единицу времени, wj - страховой запас для товара i-й номенклатуры.Минимизация критерия (37) при ограничениях (38) осуществляется по вектору z, при этом на каждом шаге пересчитываются иа(к). Найденное значение вектора z , обеспечивает минимальные издержки на временном интервале от к до k+T и, в результате выполнения ограничений (36), обеспечивает также высокую загруженность транспортного средства в соответствии с заданным коэффициентом Kr. По найденному вектору z определяется объем поставок иа(к+Т+1), и далее, по аналогии, решается задача минимизация критерия /изд (z, k +1), при ограничениях (38)(V k е [k +1, k + T +1]) и определяется новый вектор z , который используется для определения поставок в момент времени к+Т+2, и так далее.6. Моделирование адаптивной системы управления запасамиМоделирование системы управления запасами выполнено для двухноменкла-турного склада с использованием пакета Matlab 7.5.Исходные данные при моделировании были следующими:(1 - к 0 А (1(f)(45Л (040) = {25), s,(0)-Г0), st(0).(0foi J (00,35 у 0 (01,4 J 6° ^ oio yF (1,3 /h = 1, T = 50, Gmax = 150, Kr = 0,8, Ci = 1,5, С2 = 2,5, wi = 6, wi = 12,где kl =0,001, k2 =0,0005 - коэффициенты потерь. В модели спроса использовались следующие истинные значения параметров Q{ = 2,6, 92 = 2,5. На рис. 1 - 6 приведены результаты, полученные при оптимизации критерия (37) при ограничениях (38).На рис. 1 показано изменение количества каждого вида товара x1, x2. Как видно из рис. 1, количество товара постоянно меняется за счет изменения спроса и пополнения запасов.k 0 Рис. 1На рис. 2 приведены оценки параметров модели спроса (сплошная линия), пунктирными линиями обозначены истинные значения параметров.На рис. 4 приведены графики изменения спроса (сплошная линия) и оценки прогноза спроса на один такт (пунктирная линия).4k 0 Рис. 4На рис. 5 показаны реализации поставок товаров на склад.«8,1«8,220 40k 0 Рис. 52040_zНа рис. 6 показаны в виде столбиковой диаграммы вес каждой поставки.G(%) I1'GmaxGmin10020 40 Рис. 6Как видно из рисунка, в результате выполнения ограничений вида (36) в течение 50 тактов реализовано 9 поставок, при этом вес каждой поставки находится в пределах между Gmin и Gmax, обеспечивая коэффициент использования грузоподъемности транспортного средства не хуже заданного Kr. При этом был определенвектор z , минимизирующий критерий (37) с ограничениями (38). Оптимальный вектор

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