• язык CPLEX
  • язык GMPL

  • линейное программирование

    линейное программирование - это нахождение минимакса для функции f:

           f(x₁,x₂,..,xₙ)  =   Σ ck * xk
    в ограниченной области G пространства размерности n

    ограничения:

            a₁₁*x₁ + a₁₂*x₂ + .. ≤ b₁
            a₂₁*x₁ + a₂₂*x₂ + .. ≤ b₂
            ..

    набор неотрицательных переменных:

             x₁ ≥ 0
             x₂ ≥ 0
             ..  

    в канонической форме

                  c⁺ * x
                  A * x = b
                  x ≥ 0
    
    где           x     - вектор, который должен быть найден,
                  c, b  - векторы заданных коэффициентов,
                  A     - матрица заданных коэффициентов
    

    неравенства A * x ≤ b и x ≥ 0 задают ограничения, с которыми надо оптимизировать функцию f

    впервые в терминах линейного программирования в 1939г. сформулировал задачу Леонид Канторович

    в 40-х годах экономист Tjalling Koopmans тоже сформулировал некоторые классические задачи экономики в терминах линейного программирования

    в 1975г. Канторович и Koopmans получили премию шведского банка за работы в области экономики

    в случае, если все неизвестные должны быть целыми числами задача называется integer programming (IP) или integer linear programming (ILP). такие задачи считаются NP-hard

    если лишь часть переменных должна быть целыми числами, то задача называется mixed integer programming (MIP). она, как правило, тоже NP-hard

    особый случай - когда переменные могут принимать лишь два значения (бинарные переменные). такие переменные довольно часто используются при решении задач с прерывными условиями (discontinue conditions)

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

    правая часть уравнений (b-вектор), должны быть неизменяемым вектором:

            A * x ≤ B
    равносильно:
            A * x - B ≤ 0

    если исследуемое уравнение умножить на -1, то задача максимизации станет задачей минимизации (и наоборот):

            minimize      3 * x1  -  x2
    может также быть записана как:
            maximize     -3 * x1  +  x2

    для решения задач линейного программирования используют софт-решатель glsol


    glpsol


    язык CPLEX

    ключ командной строки : glpsol --cpxlp

    строка состоит из 232 символов

    комментарии

    начинаются с символа \ и идут до конца строки

    keywords

        minimize          minimum               min
        maximize          maximum               max
        subject to        such that             s.t.           st.             st
        bounds            bound
        general           generals              gen
        integer           integers              int
        binary            binaries              bin
        infinity          inf
        free
        end    

    keywords нечувствительны к регистру

    порядок секций в LP файле

    constraint definitions

    могут быть :

    by default all variables are non-negative, i.e. have zero lower bound and no upper bound. therefore definitions of these default bounds can be omitted in the bounds section

    the general , integer , and binary sections are intended to define some variables as integer or binary. all these sections are optional and needed only in case of MIP problems. if they are specified, they must follow the bounds section or, if the latter is omitted, the constraints section


    пример определения пропорций в популяции

    группа из ста детей, половину которых случайным образом (буквально кидая монетку для каждого), но "обязательно к исполнению", отобрали в спецпрограмму обучения, а другая половина шла в обычной программе

    по факту, получили, что "показатель успеха" (мы точно знаем как определить этот фактор) для прошедших спецпрограмму 75%, а для прошедших обычную программу 45%

    обозначим:
    c00 лузер, без разницы какую программу он проходил
    c01 добился успеха после спецпрограммы, но был бы лузером после обычной программы
    c10 добился успеха после обычной программы, но был бы лузером после спецпрограммы
    c11 виннер, без разницы в какой программе учился

    нам необходимо определить пропорции в популяции, приводящие к максимальному успеху в целом по популяции - при заданных условиях успешности

    код для решателя:

    maximize 
       success : c10 + c01
    
    s.t.
       b1 : c01 + c11 = 0.75 
       \ это мы имеем по результатам замеров. тогда комплиментарно мы имеем также и
       b2 : c00 + c10 = 0.25 
    
       b3 : c10 + c11 = 0.45 
       \ это мы имеем по результатам замеров. тогда комплиментарно мы имеем также и
       b4 : c00 + c01 = 0.55
    
    bounds 
       0 <= c00 <= 1
       0 <= c01 <= 1
       0 <= c10 <= 1
       0 <= c11 <= 1
    
    \ glpsol --cpxlp golos.lp -o output.txt
    

    запускаемся и смотрим файл результатов:

    Problem:    
    Rows:       4
    Columns:    4
    Non-zeros:  8
    Status:     OPTIMAL
    Objective:  success = 0.8 (MAXimum)
    
       No.   Row name   St   Activity     Lower bound   Upper bound    Marginal
    ------ ------------ -- ------------- ------------- ------------- -------------
         1 b1           NS          0.75          0.75             =             1 
         2 b2           NS          0.25          0.25             =             2 
         3 b3           NS          0.45          0.45             =            -1 
         4 b4           B           0.55          0.55             = 
    
       No. Column name  St   Activity     Lower bound   Upper bound    Marginal
    ------ ------------ -- ------------- ------------- ------------- -------------
         1 c10          B           0.25             0             1 
         2 c01          B           0.55             0             1 
         3 c11          B            0.2             0             1 
         4 c00          NL             0             0             1            -2 
    
    Karush-Kuhn-Tucker optimality conditions:
    
    KKT.PE: max.abs.err = 0.00e+00 on row 0
            max.rel.err = 0.00e+00 on row 0
            High quality
    
    KKT.PB: max.abs.err = 0.00e+00 on row 0
            max.rel.err = 0.00e+00 on row 0
            High quality
    
    KKT.DE: max.abs.err = 0.00e+00 on column 0
            max.rel.err = 0.00e+00 on column 0
            High quality
    
    KKT.DB: max.abs.err = 0.00e+00 on row 0
            max.rel.err = 0.00e+00 on row 0
            High quality
    
    End of output 


    пример с рационом

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

        c⁺ * x , xi > 0 , ci > 0 
        A * x = b  , aij > 0 , bi > 0 

    \ жиры белки углеводы в рационе при ограничении стоимости пайка
    \   c    цена продукта          c1 = 1.1    c2 = 1.3 
    \   b    жиры ; белки ; углеводы в сумме продуктов
    \   A    жиры, белки, углеводы - в каждом продукте
    \   необходимо найти вектор x  - суммарный продукт в количестве x1 x2 
    
    min
      f : 1.1 x1 + 1.3 x2 
      
    s.t.
       b1 : 1 x1 + 3 x2 > 20 
       b2 : 2 x1 + 2 x2 > 30
       b3 : 3 x1 + 2 x2 > 40
    
    bounds
      x1 > 0
      x2 > 0
    

    Problem:    
    Rows:       3
    Columns:    3
    Non-zeros:  9
    Status:     OPTIMAL
    Objective:  f = 22 (MINimum)
    
       No.   Row name   St   Activity     Lower bound   Upper bound    Marginal
    ------ ------------ -- ------------- ------------- ------------- -------------
         1 b1           B             40            30               
         2 b2           NL            40            40                         0.3 
         3 b3           NL            50            50                         0.2 
    
       No. Column name  St   Activity     Lower bound   Upper bound    Marginal
    ------ ------------ -- ------------- ------------- ------------- -------------
         1 x1           NL             0             0                        0.02 
         2 x2           B             10             0               
         3 x3           B             10             0               
    
    Karush-Kuhn-Tucker optimality conditions:
    
    KKT.PE: max.abs.err = 0.00e+00 on row 0
            max.rel.err = 0.00e+00 on row 0
            High quality
    
    KKT.PB: max.abs.err = 0.00e+00 on row 0
            max.rel.err = 0.00e+00 on row 0
            High quality
    
    KKT.DE: max.abs.err = 0.00e+00 on column 0
            max.rel.err = 0.00e+00 on column 0
            High quality
    
    KKT.DB: max.abs.err = 0.00e+00 on row 0
            max.rel.err = 0.00e+00 on row 0
            High quality
    
    End of output


    пример с двумя упаковочными фабриками

    компания по производству тары производит три вида продукции: коробки, мешки и пакеты

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

    у компании две фабрики: A и Б; каждая фабрика может работать не более, чем 24 дня в месяц, но и не менее чем 12 дней (КЗОТ, профсоюзы, все такое)

    ежедневная производительность фабрики A составляет 2 коробки, 4 мешка и 6 пакетов
    ежедневная производительность фабрики Б составляет 1 коробка, 3 мешка и 5 пакетов

    полные ежедневные затраты A составляют $100
    полные ежедневные затраты Б составляют $75

    вы - старший экономист этой компании. сколько дней в месяц должна отработать каждая фабрика, чтобы, выполнив обязательства, минимизировать расходы?

    переменные: x = дни работы фабрики A, y = дни работы фабрики Б
    целевая функция: f = 100 * x + 75 * y
    ограничения: обязательства и возможности

                    2 * x  +  1 * y ≥  50
                    4 * x  +  3 * y ≥ 120
                    6 * x  +  5 * y ≤ 250
    
                    12 ≤ x ≤ 24
                    12 ≤ y ≤ 24

    кодим

    \ переменные:       x = дни работы фабрики A,         y = дни работы фабрики Б 
    \ целевая функция:  f = 100 * x + 75 * y 
    \ ограничения:      обязательства (b1, b2) и возможности (b3)
    \                   дни работы фабрик
    
    Minimize
    f : 100 x + 75 y 
    
    Subject To 
    b1 : 2 x  +  1 y =   50
    b2 : 4 x  +  3 y =  120
    b3 : 6 x  +  5 y <= 250
    
    Bounds 
    12 <= x <= 24
    12 <= y <= 24
    
    Integer 
    x y
    

    запускаемся glpsol --cpxlp factory.lp -o output.txt и смотрим файл output.txt:

    Problem:    
    Rows:       3
    Columns:    2 (2 integer, 0 binary)
    Non-zeros:  6
    Status:     INTEGER OPTIMAL
    Objective:  f = 3000 (MINimum)
    
       No.   Row name        Activity     Lower bound   Upper bound
    ------ ------------    ------------- ------------- -------------
         1 b1                         50            50             = 
         2 b2                        120           120             = 
         3 b3                        190                         250 
    
       No. Column name       Activity     Lower bound   Upper bound
    ------ ------------    ------------- ------------- -------------
         1 x            *             15            12            24 
         2 y            *             20            12            24 
    
    Integer feasibility conditions:
    
    KKT.PE: max.abs.err = 0.00e+00 on row 0
            max.rel.err = 0.00e+00 on row 0
            High quality
    
    KKT.PB: max.abs.err = 0.00e+00 on row 0
            max.rel.err = 0.00e+00 on row 0
            High quality
    
    End of output 


    пример с тремя упаковочными фабриками

    \ компания должна произвести за следующий календарный месяц : 50 коробок и 120 мешков
    \ компания не может произвести за месяц более : 250 пакетов
    \ переменные:       дни работы: x = фабрика A , y = фабрика Б , z = фабрика В  
    \ целевая функция:  f = 100 * x + 75 * y   +  70 z
    \ ограничения:      обязательства (b1, b2) и возможности (b3)
    \ ежедневная производительность фабрики A : 2 коробки, 4 мешка и 6 пакетов
    \                                       Б : 1 коробка, 3 мешка и 5 пакетов
    \                                       В : 3 коробки, 2 мешка и 4 пакета
    \  рабочие дни фабрик                   12 - 24
    
    Minimize
    f : 100 x + 75 y + 70 z 
    
    Subject To 
    b1 : 2 x  +  1 y  +  3 z  =   50
    b2 : 4 x  +  3 y  +  2 z  =  120
    b3 : 6 x  +  5 y  +  4 z  <= 250
    
    Bounds 
    12 <= x <= 24
    12 <= y <= 24
    12 <= z <= 24
    
    Integer 
    x y z
    

    запускаемся glpsol --cpxlp factory3.lp -o output.txt и смотрим файл результатов:

    Problem:    
    Rows:       3
    Columns:    3 (3 integer, 0 binary)
    Non-zeros:  9
    Status:     INTEGER EMPTY
    Objective:  f = 0 (MINimum)
    
       No.   Row name        Activity     Lower bound   Upper bound
    ------ ------------    ------------- ------------- -------------
         1 b1                          0            50             = 
         2 b2                          0           120             = 
         3 b3                          0                         250 
    
       No. Column name       Activity     Lower bound   Upper bound
    ------ ------------    ------------- ------------- -------------
         1 x            *              0            12            24 
         2 y            *              0            12            24 
         3 z            *              0            12            24 
    
    Integer feasibility conditions:
    
    KKT.PE: max.abs.err = 0.00e+00 on row 0
            max.rel.err = 0.00e+00 on row 0
            High quality
    
    KKT.PB: max.abs.err = 1.20e+02 on row 2
            max.rel.err = 9.92e-01 on row 2
           SOLUTION IS INFEASIBLE
    
    End of output

    вот так-то, малята... решения в целых числах - нет


    язык GMPL

    ключ коммандной строки : glpsol --math

    комментарии

    начинаются с символа # и идут до конца строки

    модель

    модель задачи линейного программирования записывается на языке GMPL с использованием следующих видов объектов:

    каждый объект снабжается уникальным именем, используемым для ссылок на этот объект при построении модели

    объекты могут представлять собой многомерные массивы, адресуемые с помощью индексных множеств. для обращения к одному элементу такого объекта используется нотация object[index], где object - это имя объекта, а index - это перечисление значений индексных множеств. например, если demand является двумерным параметром, определенным на множестве IxJ, то конкретный элемент параметра demand будет адресоваться как demand[i,j], где i ∈ I и j ∈ J

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

    описание модели

    все предложения GMPL можно разделить на две категории: декларативные и функциональные

    каждое предложение должно заканчиваться точкой с запятой ;

    множество (set)

    предложение set определяет множество. если не задано пространство индексов, то set определяет простое множество, в противном случае - массив множеств по количеству элементов пространства индексов

    set name alias domain, attrib, …, attrib;
    name - символическое имя определяемого множества
    alias - необязательная строка, задающая псевдоним множества
    domain - выражение, задающее пространство индексов объекта
    attrib, …, attrib - необязательные атрибуты множества:

    примеры
    set nodes;
    это предложение определяет множество вершин, не указывая, какие элементы этому множеству принадлежат. подобная запись имеет смысл, если предполагается, что элементы множества nodes будут заданы в секции данных
    set arcs within nodes cross nodes;
    данное предложение определяет множество дуг, указывая, с помощью атрибута within, что это множество должно быть подмножеством декартова произведения множества nodes на само себя. значения элементов множества также должны быть заданы в секции данных, а атрибут within будет использован для контроля правильности заполнения этого множества
    set DaysOfWeek := {Mon, Tue, Wed, Thu, Fri, Sat, Sun};
    в этом предложении определяется множество дней недели и сразу задаются элементы множества

    параметры (param)

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

    param name alias domain, attrib, …, attrib;
    name - символическое имя определяемого параметра;
    alias - необязательная строка, задающая псевдоним параметра;
    domain - выражение, задающее пространство индексов параметра;
    attrib, …, attrib - необязательные атрибуты параметра:

    примеры
    param N;
    определение числового параметра N, который может принимать любые значения. значение параметра не определено; предполагается, что оно будет задано в секции данных
    param N integer;
    определение параметра N, значениями которого могут быть только целые числа
    param price{e in arcs};
    при определении параметра price используется пространство индексов (задаваемое множеством arcs), следовательно, price задает массив параметров - по одному числовому значению для каждого элемента множества arcs
    param demand{i in 1..N};
    еще один пример параметра с пространством индексов. в данном случае, в роли этого пространства выступает множество целых чисел от 1 до N. обращение к отдельным значениям параметра возможно с помощью явного указания индекса, например demand[3]

    переменные (var)

    предложение var определяет переменную задачи или целый массив переменных (в случае, если задано пространство индексов). значения переменные получают в ходе решения задачи оптимизации

    var name alias domain, attrib, …, attrib;
    name - символическое имя определяемой переменной
    alias - необязательная строка, задающая псевдоним переменной
    domain –выражение, задающее пространство индексов переменной
    attrib, …, attrib - необязательные атрибуты переменной:

    примеры
    var x >= 0;
    определение неотрицательной переменной
    var flow{a in arcs} >= 0;
    при определении переменной flow используется пространство индексов, следовательно, это предложение определяет целый массив переменных - по одной для каждого элемента множества arcs, причем каждая из этих переменных должна быть неотрицательной
    var take{o in objects} binary;
    задает массив бинарных переменных по одной для каждого элемента множества objects

    ограничения (subject to)

    предложение subject to определяет ограничение или целую группу ограничений (в случае, если задано пространство индексов). ключевое слово subject to может быть записано сокращенно как subj to, s.t. или даже пропущено

    subject to name alias domain : expression, = expression ;
    subject to name alias domain : expression, ≤ expression ;
    subject to name alias domain : expression, >= expression ;
    subject to name alias domain : expression, ≤ expression, ≤ expression;
    subject to name alias domain : expression, >= expression, >= expression;

    name - символическое имя определяемого ограничения;
    alias - необязательная строка, задающая псевдоним ограничения;
    domain –выражение, задающее пространство индексов ограничения;
    expression - линейное выражение, используемое для вычисления компонента ограничения

    примеры
    s.t. r: x[1] + x[2] ≤ 17;
    данное ограничение требует, чтобы сумма двух переменных не превосходила 17. из самой записи не следует, что x - это переменные, однако это предполагается. следует заметить, что подобное ограничение является образцом плохого стиля, поскольку константа 17 присутствует непосредственно в определении модели - правильней было бы определить соответствующий параметр и использовать его в этом ограничении
    subject to weight: sum{o in objects} w[o] * take[o] ≤ W_max;
    ограничение, которое могло бы быть использовано в задаче об упаковке рюкзака. при такой интерпретации weight - это имя ограничения, выражение sum используется для получения суммарного веса набора объектов, для которых переменная take принимает значение 1, W_max - параметр, задающий вместительность рюкзака, а w - параметр-массив (индексируемый элементами множества objects), в котором содержатся веса объектов
    subject to ration{p in products}: sum{c in crops} x[c] * r[p, c] ≤ R[p];
    это предложение задает сразу семейство однотипных ограничений - по одному для каждого элемента множества products

    целевая функция (maximize/minimize)

    предложение minimize/maximize определяет целевую функцию. при добавлении в модель нескольких целевых функций, использоваться в процессе решения задачи будет только первая из них

    maximize name alias : expression ;
    minimize name alias : expression ;

    name - символическое имя целевой функции;
    alias - необязательная строка, задающая псевдоним целевой функции;
    expression - линейное выражение, используемое для вычисления значения целевой функции

    примеры
    minimize total_cost: sum{w in warehouses, c in consumers} c[w, c] * x[w, c] ;
    пример целевой функции, которая могла бы быть использована в транспортной задаче, если множество пунктов хранения обозначено warehouses, множество потребителей - consumers, цена перевозки единицы груза обозначена параметром c, а объем перевозок - переменной x
    maximize total_cost: sum{o in objects} price[o] * take[o] ;
    пример целевой функции, которая могла бы быть использована в задаче о рюкзаке

    решение задачи (solve)

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

    вывод найденного решения (display, printf, table)

    для вывода в GMPL предусмотрено несколько функциональных предложений. самое простое из них - display. для display может быть задано пространство индексов, а выводимые значения разделяются запятой:

            display : ‘x= ‘, x, ‘y = ‘, y ;
            display{o in objects} : o, take[o] ; 

    с помощью display можно выводить не только значения выражений, но и объекты модели (множества, ограничения, целевую функцию), однако вывод всегда производится в стандартный поток (как правило, консоль) и возможности форматирования существенно ограничены

    большую гибкость предоставляет предложение printf, которое позволяет осуществлять форматированный вывод, используя синтаксис описания формата одноименной функции языка C. кроме того, printf позволяет направлять вывод в указанный файл (как с перезаписью (“>”), так и с добавлением в конец файла (“>>”)):

            printf : "Hello, world!\n" ; 
            printf{i in 1..4} : “x[%d] = %.3f\n”, i, x[i] > “result.txt” ;

    предусмотрен и высокоуровневый интерфейс для ввода-вывода табличных данных с помощью предложения table:

       table name alias IN driver arg … arg: set <- [fld, …, fld] , par ~ fld, …, par ~ fld ;
       table name alias domain OUT driver arg … arg : expr ~ fld, … , expr ~ fld ;
    name –имя таблицы
    alias - необязательная строка, задающая псевдоним таблицы
    domain - выражение, задающее пространство индексов таблицы;
    IN - означает, что данные должны быть прочитаны из таблицы
    OUT - означает, что данные должны быть записаны в таблицу
    driver - идентификатор драйвера, предоставляющего доступ к таблице (поддерживаются файлы формата CSV, базы данных ODBC, xBASE, MySQL)
    arg - необязательные параметры, которые должны быть переданы драйверу (например, имя файла для CSV)
    set - имя контрольного множества
    fld - имя поля
    par - имя параметра модели (пространство индексов этого параметра должно содержать контрольное множество)
    expr - числовое или символьное выражение

    описание данных

    описание данных расположено либо в отдельной секции файла модели, начинающейся предложением data ; и заканчивающимся предложением end ; , либо в отдельном файле (этот файл не обязательно начинать с предложения data ; )

    описание данных должно задавать элементы всех множеств, определенных в секции описания модели, и значения всех указанных там параметров

    задание элементов множества
    set name, record, …, record ;
    set name [symbol, …, symbol] , record, …, record ;

    name - имя множества;
    symbol, …, symbol - индексы, задающие конкретный элемент, если у множества соответствующего имени name, определено пространство индексов;
    record, …, record - записи с данными:

    примеры
    set nodes := LED VKO KGD TJM ;
    в этом примере элементы множества nodes задаются в простой форме. элементами этого множества являются символьные константы (являющиеся международными кодами аэропортов)
    set arcs := (LED, VKO) (VKO, LED) (VKO, KGD) (VKO, TJM) ;
    если множество arcs было определено как подмножество декартова произведения множества nodes самого на себя, то элементами этого множества должны являться пары, что и имеет место в задании arcs выше
              set arcs :     LED VKO KGD TJM :=
                         LED  -   +   -   -
                         VKO  +   -   +   +
                         KGD  -   -   -   -
                         TJM  -   -   -   - ;
    А это пример задания элементов множества пар в матричной форме. содержимое множества arcs будет таким же, как в предыдущем примере

    задание значений параметров

    param name, record, …, record ;
    param name default value, record, …, record ;

    name - имя параметра;
    value - необязательное значение по умолчанию для параметра;
    record, …, record - записи с данными:

    примеры
    : param N := 5 ;
    param demand := 1 8 2 14 3 3 4 9 5 11 ;
    параметр N был определен как скалярный параметр, задание значения для него имеет очевидный синтаксис. параметр demand определен как параметр-массив, поскольку для него было задано индексное пространство 1..N. и поскольку N присвоено значение 5, то demand должен состоять из пяти значений, которые и перечисляются с указанием соответствующих индексов. то есть, demand[1] присвоено значение 8, а demand[5] - значение 11
    param r : 1 2 3 := 1 0.5 0.3 0.12 2 0.2 0.1 0.22 ;
    это пример использования записи tabular-data для задания параметра с двумерными индексами. при таком задании r[2,1] будет присвоено значение 0.2
    param trans_cost :=
    [*,*,wool] :
         1     2     3 :=
    1   0.3   0.2   0.11
    2   0.1   0.1   0.32
    [*,*,steel]:
         1     2     3 :=
    1   0.7   0.4   0.45
    2   0.2   0.2   0.64 ; 
    такая конструкция потребуется, если необходимо задать значения для параметра, пространство индексов которого образовано кортежами из трех элементов. например, стоимость перевозки единицы товара зависит не только от того, откуда и куда товар нужно перевезти, но и от самого товара. в данном примере, товаров рассматривается два: шерсть (wool) и сталь (steel). конструкция [*,*,wool] называется спецификация среза и говорит о том, что следующий за ней блок данных следует воспринимать как задание двумерного параметра (по количеству символов “*”) описывающего стоимости перевозки шерсти (wool). то есть, задание значений многомерного параметра происходит по срезам, размерность каждого из которых не превышает 2

    пример : первая задача Беллмана

    пример решения первой задачи Беллмана при использовании марковских цепей. все в одном файле. если ключ -o не указан в командной строке, то файл с решением создан не будет и единственный способ увидеть решение - это использовать display, printf или table

    #                glpsol -m markov.mod       
    
    set States ;
    set Actions ; 
    
    var value {s in States} ; 
    
    param transition_probability {s in States , a in Actions , s2 in States} ; 
    param reward {s in States , a in Actions} ;
    param discount_factor ;
    
    minimize
    total: sum {s in States} value[s] ;
    
    s.t.
      bellman {s in States , a in Actions}: 
         value[s] >= 
           reward[s,a] + 
           sum {s2 in States} discount_factor * transition_probability[s , a , s2] * value[s2] ;
    
    solve ;
    
    display : total ;
    
    end ;
            
    data ; 
    
    set States := good deteriorating broken ;
    set Actions := maintain ignore ;
    
    param transition_probability:= 
      [* , * , good]: 
                    maintain ignore := 
        good            1      .5 
        deteriorating  .9       0 
        broken         .2       0 
      [* , * , deteriorating]: 
                   maintain ignore := 
        good            0      .5 
        deteriorating  .1      .5 
        broken          0       0 
      [* , * , broken]: 
                   maintain ignore := 
        good           0        0 
        deteriorating  0       .5 
        broken        .8        1 ;
    
    param reward: 
                   maintain ignore := 
        good            1      2 
        deteriorating   1      2 
        broken         -1      0 ; 
        
    param discount_factor := .8 ; 
    
    end ;
    


    пример трех цехов на одной фабрике

    фабрика производит три вида продукции: P1, P2 и P3. известна цена на продукцию и спрос на каждый из видов продукции

    продукт цена спрос
    P1 1200 35
    P2 2500 25
    P3 1400 30

    на фабрике есть три цеха: W1, W2 и W3. для производства продукции P1 необходимы только технологические операции, производимые цехом W1, для P2 - W1 и W3, для производства P3 необходима обработка во всех трех цехах W1, W2 и W3. причем, если в цехах W1 и W2 продукция разных видов обрабатывается одинаково, и известна общая производительность этих цехов в единицах обработанной продукции, то в цехе W3 известны как временные затраты на обработку каждого из видов продукции, так и общий фонд рабочего времени сотрудников цеха

    W1, шт. W2, шт. W3, часов (P2/P3/фонд)
    40 20 8/2/80

    из всех видов материалов, используемых при производстве продукции, ограниченным является только один, поставки его в неделю и потребности для каждого из видов продукции приведены в таблице

    поставки, кг на ед. P1, кг на ед. P2, кг на ед. P3, кг
    50 0,8 0,6 0,7

    необходимо составить производственный план, максимизирующий выручку от реализации продукции

    файл модели, factory1.mod :

    set P := 1..3 ; 
    set W := 1..3 ; 
    set W_Mass within W := {1, 2} ; 
    set W_Manual within W := W diff W_Mass ;
    set Ops within P cross W := {(1, 1), (2, 1), (2, 3), (3, 1), (3, 2), (3, 3)} ;
    set W_P {w in W} := setof {(i, w) in Ops} i ;
    
    param price {p in P} ;                      # цена по видам продукции, руб. 
    param demand {p in P} ;                     # спрос в неделю по видам продукции, шт. 
    param cost {p in P} ;                       # расход материала на производство продукции, кг 
    param prod {w in W_Mass} ;                  # производительность цехов, шт.   
    param prod_m_cost {w in W_Manual, p in P} ; # затраты времени на единицу продукции, часы 
    param prod_m {w in W_Manual} ;              # фонд рабочего времени, часы 
    param materials ;                           # поставки материалов в неделю, кг
    
    var x {p in P} >= 0 ;
    
    # максимальное количество изделий, которые могут пройти машинную обработку :
    s.t. 
    mass {w in W_Mass} : sum {p in W_P[w]} x[p] <= prod[w] ;
    
    # для цеха с ручной обработкой должно выполняться требование по рабочему времени : 
    s.t. 
    manual {w in W_Manual} : sum {p in W_P[w]} x[p] * prod_m_cost[w, p] <= prod_m[w] ;
    
    # производство каждого из видов продукции не должно превосходить спрос : 
    s.t. 
    dem {p in P} : x[p] <= demand[p] ;
    
    # затраты материала на производство не должны превышать поставок : 
    s.t.
    mat : sum {p in P} x[p] <= materials ;
    
    maximize 
    profit : sum {p in P} price[p] * x[p] ;
    
    solve ; 
    
    end ;
    
    # > glpsol -m factory1.mod -d factory1.dat -o factory1.sol 

    файл данных, factory1.dat:

    data ; 
    
    # цена по видам продукции, руб. 
    param price := 1 1200 2 2500 3 1400 ; 
    
    # спрос по видам продукции, шт. 
    param demand := 1 35 2 25 3 30 ; 
    
    # расход материала на производство продукции, кг 
    param cost := 1 0.8 2 0.6 3 0.7 ; 
    
    # производительность цехов, шт. 
    param prod := 1 40 2 20 ; 
    
    # затраты времени на единицу продукции, часы 
    param prod_m_cost default 0 := [3, 2] 8 [3, 3] 2 ; 
    
    # фонд рабочего времени, часы 
    param prod_m := 3 80 ; 
    
    # поставки материалов, кг 
    param materials := 50 ; 
    
    end ;
    

    запускаемся : glpsol -m factory1.mod -d factory1.dat -o factory1.sol

    файл результатов, factory1.sol:

    Problem:    factory1
    Rows:       8
    Columns:    3
    Non-zeros:  15
    Status:     OPTIMAL
    Objective:  profit = 61000 (MAXimum)
    
       No.   Row name   St   Activity     Lower bound   Upper bound    Marginal
    ------ ------------ -- ------------- ------------- ------------- -------------
         1 mass[1]      NU            40                          40          1200 
         2 mass[2]      B              0                          20 
         3 manual[3]    NU            80                          80         162.5 
         4 dem[1]       B             30                          35 
         5 dem[2]       B             10                          25 
         6 dem[3]       B              0                          30 
         7 mat          B             40                          50 
         8 profit       B          61000                             
    
       No. Column name  St   Activity     Lower bound   Upper bound    Marginal
    ------ ------------ -- ------------- ------------- ------------- -------------
         1 x[1]         B             30             0               
         2 x[2]         B             10             0               
         3 x[3]         NL             0             0                        -125 
    
    Karush-Kuhn-Tucker optimality conditions:
    
    KKT.PE: max.abs.err = 0.00e+00 on row 0
            max.rel.err = 0.00e+00 on row 0
            High quality
    
    KKT.PB: max.abs.err = 0.00e+00 on row 0
            max.rel.err = 0.00e+00 on row 0
            High quality
    
    KKT.DE: max.abs.err = 0.00e+00 on column 0
            max.rel.err = 0.00e+00 on column 0
            High quality
    
    KKT.DB: max.abs.err = 0.00e+00 on row 0
            max.rel.err = 0.00e+00 on row 0
            High quality
    
    End of output 

    оптимальное значение целевой функции (строка profit - по имени целевой функции, введенному в файле модели) 61000 достигается при производстве 30 штук продукции первого вида (значение из столбца Activity напротив переменной, соответствующей плану выпуска первой продукции - x[1]) и 10 штук продукции второго вида (аналогично, для переменной x[2]). Переменная x[3] помечена статусом NL (небазисная переменная, значение которой находится на нижней границе); действительно, в решении эта переменная принимает значение 0 (столбец Activity), а это является нижней границей области допустимых значений неотрицательной переменной. поскольку данная переменная является небазисной, она имеет ненулевую приведенную цену -125 (столбец Marginal). то есть, увеличение значения переменной x[3] на единицу повлечет за собой снижение значения целевой функции на 125 единиц

    судя по значениям столбца St первой таблицы, из всех ограничений задачи активными являются два: mass[1] (ограничение производительности первого цеха) и manual[3] (ограничение производительности третьего цеха). статус NU (неравенство с активным ограничением сверху) согласуется с тем, что значение левой части ограничения (столбец Activity) для данных ограничений совпадает с верхней границей (столбец Upper bound). столбец Marginal содержит теневые цены для этих ограничений (1200 и 162,5 соответственно)

    анализ чувствительности

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

    для создания такого отчета при вызове glpsol следует указать параметр --ranges, задав имя файла, куда отчет должен быть помещен. например:

          $> glpsol -m factory1.mod -d factory1.dat --ranges factory1_sens.sol -o factory1.sol    

    анализ чувствительности активных границ проводится только для активных ограничений. для каждого такого ограничения в ходе исследования активная граница (правая часть ограничения) изменяется как в сторону уменьшения, так и в сторону увеличения. каждому ограничению в отчете соответствует две строки: первая характеризует результат, полученный при уменьшении активной границы, вторая - при увеличении. изменение активной границы ведет к изменению значений переменных и продолжается до тех пор, пока какая-либо из базисных переменных не достигает границы области допустимых значений. соответствующее значение активной границы называется в GLPK break point и по сути является границей интервала осуществимости решения по отношению к исследуемому ограничению. в отчете по анализу чувствительности диапазон исследуемых значений границы приводится в столбце Activity range. соответствующие значения целевой функции указаны в столбце Obj value at break point, а символические имена ограничивающих базисных переменных - в столбце Limiting variable. если активную границу можно увеличивать или уменьшать бесконечно, то в столбце Activity range будет указано +Inf или –Inf соответственно, а столбец Limiting variable будет пуст

    файл factory1_sens.sol:

    GLPK 4.65 - SENSITIVITY ANALYSIS REPORT                    Page   1
    
    Problem:    factory1
    Objective:  profit = 61000 (MAXimum)
    
     No. Row name St Activity Slack Lower bound Activity Obj coef Obj value at Limiting
       Marginal Upper bound range range break point variable
    ---------- ------------- ------------- ------------- ------------- ------------- ------------
     1 mass[1] NU 40.00000 .  -Inf 10.00000 -1200.00000 25000.00000 x[1] 
       1200.00000 40.00000 45.00000 +Inf 67000.00000 dem[1]
    
     2 mass[2] BS .  20.00000 -Inf .  -Inf 61000.00000
       .  20.00000 30.00000 125.00000 61000.00000 x[3]
    
     3 manual[3] NU 80.00000 .  -Inf 40.00000 -162.50000 54500.00000 dem[1]
       162.50000 80.00000 200.00000 +Inf 80500.00000 dem[2]
    
     4 dem[1] BS 30.00000 5.00000 -Inf 15.00000 -166.66667 56000.00000 x[3]
       .  35.00000 40.00000 1300.00000 100000.00000 manual[3]
    
     5 dem[2] BS 10.00000 15.00000 -Inf 5.00000 -500.00000 56000.00000 x[3]
       .  25.00000 10.00000 +Inf +Inf
    
     6 dem[3] BS .  30.00000 -Inf .  -Inf 61000.00000
       .  30.00000 20.00000 125.00000 61000.00000 x[3]
    
     7 mat BS 40.00000 10.00000 -Inf 10.00000 -1200.00000 13000.00000 mass[1]
       .  50.00000 40.00000 +Inf +Inf
    
     8 profit BS 61000.00000 -61000.00000 -Inf 25000.00000 -1.00000 . mass[1]
       .  +Inf 61000.00000 +Inf +Inf
    
    GLPK 4.65 - SENSITIVITY ANALYSIS REPORT                       Page 2
    
    Problem: factory1
    Objective: profit = 61000 (MAXimum)
    
     No. Column name St Activity Obj coef Lower bound  Activity Obj coef Obj value at Limiting
           Marginal Upper bound  range  range break point variable
    ----- ------------- ------------- ------------- ------------- ------------- ------------
      1 x[1]  BS 30.00000 1200.00000  .  15.00000 1033.33333 56000.00000 x[3]
            .   +Inf  35.00000 2500.00000 100000.00000 manual[3]
    
      2 x[2]  BS 10.00000 2500.00000  .   5.00000 2000.00000 56000.00000 x[3]
            .   +Inf  10.00000  +Inf  +Inf
    
      3 x[3]  NL  .  1400.00000  .  -6.66667  -Inf 61833.33333 dem[1]
           -125.00000  +Inf  20.00000 1525.00000 58500.00000 mass[2]
    
    End of report

    например, ограничение mass[1] является активным в верхней границе (статус NU), по условию задачи равной 40 (Upper bound). соответственно, в ходе анализа чувствительности эта граница изменялась в сторону уменьшения до тех пор, пока при значении 10 переменная x[1] не приняла граничного значения (стала равной нулю), и в сторону увеличения до тех пор, пока при значении 45 не активизировалось ограничение dem[1] (спрос на продукцию первого вида). соответствующие значения целевой функции - 25000 и 67000

    анализ чувствительности коэффициентов целевой функции при небазисных переменных является достаточно простым, поскольку изменение коэффициента целевой функции при небазисной переменной ведет к эквивалентному изменению ее приведенной цены. для каждой небазисной переменной коэффициент при этой переменной в целевой функции изменяется как в сторону уменьшения, так и в сторону увеличения. каждой переменной в отчете соответствует две строки: первая описывает результат, полученный при уменьшении коэффициента, вторая - при увеличении. изменение коэффициента ведет к изменению приведенной цены анализируемой переменной и влияет на приведенные цены других переменных. базисное решение остается допустимым и оптимальным до тех пор, пока приведенная стоимость сохраняет свой знак. в отчете по анализу чувствительности диапазон значений коэффициента выводится в столбце Obj coef range. если коэффициент можно увеличивать или уменьшать бесконечно без смены базиса, то в столбце Obj coef range будет указано +Inf или –Inf соответственно

    например, переменная x[3] - небазисная переменная, принимающая в решении наименьшее допустимое значение (статус NL). коэффициент при этой переменной в исходной задаче равен 1400, а приведенная стоимость равна -125 (столбец Obj coef/Marginal). при уменьшении коэффициента при x[3] потенциальный вклад x[3] в целевую функцию становится еще меньше; что отражается и в значении -Inf в столбце Obj coef range. то есть, можно сколько угодно уменьшать данный коэффициент без изменения базиса

    уменьшение цены на третий вид продукции не приведет к пересмотру производственного плана, в частности, к включению этого вида продукции в план. при увеличении же коэффициента до 1525 состав базисных переменных изменяется и разница между исходным значением коэффициента и пороговым равна 125, то есть, объяснимо совпадает с приведенной ценой переменной x[3]

    анализ чувствительности коэффициентов целевой функции при базисных переменных - для каждой базисной переменной коэффициент при этой переменной в целевой функции изменяется как в сторону уменьшения, так и в сторону увеличения. каждой переменной в отчете соответствует две строки: первая описывает результат, полученный при уменьшении коэффициента, вторая - при увеличении. изменение коэффициента ведет к изменению приведенной цены небазисных переменных. при определенном изменении какая-то из цен становится равной нулю, а это означает, что даже малейшее дальнейшее изменение приведет к включению в базис соответствующей переменной. в отчете по анализу чувствительности диапазон значений коэффициента выводится в столбце Obj coef range. соответствующие значения целевой функции приведены в столбце Obj value at break point, символические имена ограничивающих небазисных переменных - в столбце Limiting variable, а значения исследуемой переменной, которые она принимает на концах диапазона - в столбце Activity range. если коэффициент можно увеличивать или уменьшать бесконечно без смены базиса, то в столбце Obj coef range будет указано +Inf или –Inf соответственно

    например, переменная x[1] является базисной и принимает в решении значение 30 (столбец Activity), а коэффициент в целевой функции при этой переменной равен 1200 (столбец Obj coef/Marginal). анализ показал, что переменная остается базисной при коэффициенте в целевой функции не ниже 1033.33333, при дальнейшем снижении коэффициента x[1] выйдет из базиса и ее место займет x[3]. значение переменной x[1] при граничном значении коэффициента будет равно 15 (столбец Activity range), а значение целевой функции – 56000 (столбец Obj value at break point). увеличивать коэффициент при x[1] без смены базиса можно вплоть до значения 2500. значения переменных при этом будут изменяться, значение целевой функции возрастет до 100000, но состав ненулевых переменных и активных ограничений будет неизменен. дальнейшее повышение коэффициента ограничено manual[3]

    ограничение manual[3] в решении является активным. это ограничение является ограничением нестрогого неравенства, а значит, при переходе к канонической форме в него должна быть добавлена дополнительная переменная, значение которой равно разнице между левой и правой частями исходного ограничения. при активном ограничении левая и правая части равны, а значит, дополнительная переменная не содержится в базисе. при повышении коэффициента целевой функции при x[1] выше 2500 дополнительная переменная входит в базис, то есть, левая его часть соответствующего ограничения становится строго меньше правой. это означает, что при ценах на продукцию первого вида выше 2500 план производства изменяется таким образом (в пользу продукции первого вида), что третий цех становится недогруженным