Исследование линейных несходимостей

Открыть сценарий в реальном времени

В этом примере показано, как исследовать линейные ограничения, которые приводят к невозможности решения проблемы. Для получения дополнительной информации об этих методах см. Chinneck [1] и [2].

Если линейные ограничения делают проблему неосуществимой, может понадобиться найти подмножество ограничений, которое невозможно, но удаление любого члена подмножества делает возможной оставшуюся часть подмножества. Имя такого подмножества - Неприводимое несбыточное подмножество ограничений, сокращенно IIS. И наоборот, может потребоваться найти максимально возможное подмножество ограничений. Это подмножество называется Максимально выполнимым подмножеством, сокращенно MaxFS. Эти две концепции взаимосвязаны, но не идентичны. Проблема может иметь множество различных IIS, некоторые с различной кардинальностью.

В этом примере показаны два способа поиска IIS и два способа получения возможного набора ограничений. В примере предполагается, что все заданные границы верны, что означает lb и ub аргументы не содержат ошибок.

Недостижимый пример

Создание случайной матрицы A представление линейных неравенств размера 150-на-15. Установка соответствующего вектора b к вектору с записями 10 и измените 5% от этих значений на -10.

N = 15;
rng default
A = randn([10*N,N]);
b = 10*ones(size(A,1),1);
Aeq = [];
beq = [];
b(rand(size(b)) <= 0.05) = -10;
f = ones(N,1);
lb = -f;
ub = f;

Проверьте, что problem неосуществим.

[x,fval,exitflag,output,lambda] = linprog(f,A,b,Aeq,beq,lb,ub);

No feasible solution found.

Linprog stopped because no point satisfies the constraints.

Фильтр удаления

Чтобы определить IIS, выполните следующие действия. Учитывая набор линейных ограничений с номерами от 1 до N, где все проблемные ограничения неосуществимы:

Для каждого i от 1 до N:

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

Перейти к следующему i (до значения N).

В конце этой процедуры ограничения, которые остаются в проблеме, образуют IIS.

Код MATLAB ®, реализующий эту процедуру, см. в разделеdeletionfilter вспомогательная функция в конце этого примера.

Примечание.Если для этого примера используется файл сценария в реальном времени, deletionfilter функция уже включена в конец файла. В противном случае необходимо создать эту функцию в конце файла .m или добавить ее в качестве файла по пути MATLAB. То же самое относится и к другим вспомогательным функциям, используемым ниже в этом примере.

Посмотреть эффект от deletionfilter на примере данных.

[ineqs,eqs,ncall] = deletionfilter(A,b,Aeq,beq,lb,ub);

Проблема не имеет ограничений равенства. Найти индексы для ограничений неравенства и значение b(iis).

iis = find(ineqs)

iis = 114

b(iis)

ans = -10

Только одно ограничение неравенства делает проблему неосуществимой, наряду с ограничивающими ограничениями. Ограничение:

A(iis,:)*x <= b(iis).

Почему это ограничение неосуществимо вместе с границами? Найти сумму абсолютных значений этой строки A.

disp(sum(abs(A(iis,:))))

    8.4864

Из-за границ, x вектор имеет значения от -1 до 1, и так A(iis,:)*x не может быть меньше b(iis) = –10.

Сколько их linprog звонки сделали deletionfilter выступать?

disp(ncall)

Задача имеет 150 линейных ограничений, поэтому функция называется linprog 150 раз.

Эластичный фильтр

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

Сначала разрешите каждое ограничение i нарушается на положительную сумму e(i), где ограничения равенства имеют как аддитивные, так и вычитаемые положительные значения упругости.

$\begin{array}{l} _{}_{} \\ _{}_{}_{}_{Aineqx≤bineq+eAeqx=beq+e1-e2} \end{array}$

Затем решите связанную проблему линейного программирования (LP)

$\underset{}{}_{minx,e∑ei}$

с учетом перечисленных ограничений и с $_{} ei≥0$ .

Если связанный ЛП имеет решение, удалите все ограничения, которые имеют строго положительный связанный $_{ei}$ , и запишите эти ограничения в список индексов (потенциальных членов IIS). Вернитесь к предыдущему шагу для решения нового сокращенного связанного ЛП.
Если связанный ЛП не имеет решения (неосуществимо) или не имеет строго положительного связанного $_{ei}$ , выйдите из фильтра.

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

Код MATLAB ®, реализующий этот фильтр, см. в разделеelasticfilter вспомогательная функция в конце этого примера.

Посмотреть эффект от elasticfilter на примере данных.

[ineqselastic,eqselastic,ncall] = ...
    elasticfilter(A,b,Aeq,beq,lb,ub);

Проблема не имеет ограничений равенства. Найдите индексы ограничений неравенства.

iiselastic = find(ineqselastic)

iiselastic = 5×1

     2
    60
    82
    97
   114

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

A_1 = A(ineqselastic > 0,:);
b_1 = b(ineqselastic > 0);
[dineq_iis,deq_iis,ncall2] = deletionfilter(A_1,b_1,Aeq,beq,lb,ub);
iiselasticdeletion = find(dineq_iis)

iiselasticdeletion = 5

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

disp(ncall + ncall2)

Удалить IIS в цикле

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

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

Код отслеживает исходные переменные в задаче с помощью логического вектора activeA для представления текущих ограничений (строк) A матрица и логический вектор activeAeq для представления текущих ограничений Aeq матрица. При добавлении или удалении ограничений код индексирует в A или Aeq чтобы номера строк не изменялись, даже если количество ограничений изменяется.

Выполнение этого кода возвращает idx2, вектор индексов ненулевых элементов в activeA:

idx2 = find(activeA)

Предположим, что var - логический вектор, имеющий ту же длину, что и idx2. Тогда

idx2(find(var))

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

opts = optimoptions('linprog','Display',"none");
activeA = true(size(b));
activeAeq = true(size(beq));
[~,~,exitflag] = linprog(f,A,b,Aeq,beq,lb,ub,opts);
ncl = 1;
while exitflag < 0
    [ineqselastic,eqselastic,ncall] = ...
        elasticfilter(A(activeA,:),b(activeA),Aeq(activeAeq,:),beq(activeAeq),lb,ub);
    ncl = ncl + ncall;
    idxaa = find(activeA);
    idxae = find(activeAeq);
    tmpa = idxaa(find(ineqselastic));
    tmpae = idxae(find(eqselastic));
    AA = A(tmpa,:);
    bb = b(tmpa);
    AE = Aeq(tmpae,:);
    be = beq(tmpae);
    [ineqs,eqs,ncall] = ...
        deletionfilter(AA,bb,AE,be,lb,ub);
    ncl = ncl + ncall;
    activeA(tmpa(ineqs)) = false;
    activeAeq(tmpae(eqs)) = false;
    disp(['Removed inequalities ',int2str((tmpa(ineqs))'),' and equalities ',int2str((tmpae(eqs))')])
    [~,~,exitflag] = ...
        linprog(f,A(activeA,:),b(activeA),Aeq(activeAeq,:),beq(activeAeq),lb,ub,opts);
    ncl = ncl + 1;
end

Removed inequalities 114 and equalities 
Removed inequalities 97 and equalities 
Removed inequalities 64  82 and equalities 
Removed inequalities 60 and equalities

fprintf('Number of linprog calls: %d\n',ncl)

Number of linprog calls: 28

Обратите внимание, что цикл удаляет неравенства 64 и 82 одновременно, что указывает на то, что эти два ограничения образуют IIS.

Найти MaxFS

Другой подход для получения возможного набора ограничений заключается в непосредственном поиске MaxFS. Как объясняет Чиннек [1], поиск MaxFS является NP-полной проблемой, что означает, что проблема не обязательно имеет эффективные алгоритмы для поиска MaxFS. Однако Чиннек предлагает некоторые алгоритмы, которые могут работать эффективно.

Используйте алгоритм Чиннека 7.3, чтобы найти набор ограничений обложки, который при удалении дает выполнимый набор. Алгоритм реализован в generatecover вспомогательная функция в конце этого примера.

[coversetineq,coverseteq,nlp] = generatecover(A,b,Aeq,beq,lb,ub)

coversetineq = 5×1

   114
    97
    60
    82
     2

coverseteq =

     []

nlp = 40

Удалите эти ограничения и решите проблему ЛП.

usemeineq = true(size(b));
usemeineq(coversetineq) = false; % Remove inequality constraints
usemeeq = true(size(beq));
usemeeq(coverseteq) = false; % Remove equality constraints
[xs,fvals,exitflags] = ...
    linprog(f,A(usemeineq,:),b(usemeineq),Aeq(usemeeq),beq(usemeeq),lb,ub);

Optimal solution found.

Обратите внимание, что набор обложек совпадает с набором обложек iiselastic из упругого фильтра. В общем, эластичный фильтр находит слишком большой набор крышек. Алгоритм Чиннека 7.3 начинается с результата упругого фильтра и затем сохраняет только необходимые ограничения.

Алгоритм Чиннека 7.3 принимает 40 вызовов linprog для завершения расчета MaxFS. Это число немного больше 28 вызовов, использовавшихся ранее в процессе удаления IIS в цикле.

Кроме того, обратите внимание, что неравенства, удаленные в цикле, не совсем такие же, как неравенства, удаленные алгоритмом 7.3. Цикл удаляет неравенства 114, 97, 82, 60 и 64, в то время как алгоритм 7.3 удаляет неравенства 114, 97, 82, 60 и 2. Проверьте, что неравенства 82 и 64 образуют IIS (как указано в разделе Удаление IIS в цикле), и что неравенства 82 и 2 также образуют IIS.

usemeineq = false(size(b));
usemeineq([82,64]) = true;
ineqs = deletionfilter(A(usemeineq,:),b(usemeineq),Aeq,beq,lb,ub);
disp(ineqs)

   1
   1

usemeineq = false(size(b));
usemeineq([82,2]) = true;
ineqs = deletionfilter(A(usemeineq,:),b(usemeineq),Aeq,beq,lb,ub);
disp(ineqs)

   1
   1

Ссылки

[1] Чиннек, Дж. У. Выполнимость и несходимость в оптимизации: алгоритмы и вычислительные методы. Спрингер, 2008.

[2] Чиннек, Дж. У. «Осуществимость и несходимость в оптимизации». Учебное пособие для CP-AI-OR-07, Брюсселя, Бельгии. Доступно по адресу https://www.sce.carleton.ca/faculty/chinneck/docs/CPAIOR07InfeasibilityTutorial.pdf.

Вспомогательные функции

Этот код создает deletionfilter функция помощника.

function [ineq_iis,eq_iis,ncalls] = deletionfilter(Aineq,bineq,Aeq,beq,lb,ub)
ncalls = 0;
[mi,n] = size(Aineq); % Number of variables and linear inequality constraints
f = zeros(1,n);
me = size(Aeq,1); % Number of linear equality constraints
opts = optimoptions("linprog","Algorithm","dual-simplex","Display","none");

ineq_iis = true(mi,1); % Start with all inequalities in the problem
eq_iis = true(me,1); % Start with all equalities in the problem

for i=1:mi
    ineq_iis(i) = 0; % Remove inequality i
    [~,~,exitflag] = linprog(f,Aineq(ineq_iis,:),bineq(ineq_iis),...
                                Aeq,beq,lb,ub,[],opts);
    ncalls = ncalls + 1;
    if exitflag == 1 % If now feasible
        ineq_iis(i) = 1; % Return i to the problem
    end
end
for i=1:me
    eq_iis(i) = 0; % Remove equality i
    [~,~,exitflag] = linprog(f,Aineq,bineq,...
                                Aeq(eq_iis,:),beq(eq_iis),lb,ub,[],opts);
    ncalls = ncalls + 1;
    if exitflag == 1 % If now feasible
        eq_iis(i) = 1; % Return i to the problem
    end
end
end

Этот код создает elasticfilter функция помощника.

function [ineq_iis,eq_iis,ncalls,fval0] = elasticfilter(Aineq,bineq,Aeq,beq,lb,ub)
ncalls = 0;
[mi,n] = size(Aineq); % Number of variables and linear inequality constraints
me = size(Aeq,1);
Aineq_r = [Aineq -1.0*eye(mi) zeros(mi,2*me)];
Aeq_r = [Aeq zeros(me,mi) eye(me) -1.0*eye(me)]; % Two slacks for each equality constraint
lb_r = [lb(:); zeros(mi+2*me,1)];
ub_r = [ub(:); inf(mi+2*me,1)];
ineq_slack_offset = n;
eq_pos_slack_offset = n + mi;
eq_neg_slack_offset = n + mi + me;
f = [zeros(1,n) ones(1,mi+2*me)];
opts = optimoptions("linprog","Algorithm","dual-simplex","Display","none");
tol = 1e-10;

ineq_iis = false(mi,1);
eq_iis = false(me,1);
[x,fval,exitflag] = linprog(f,Aineq_r,bineq,Aeq_r,beq,lb_r,ub_r,[],opts);
fval0 = fval;
ncalls = ncalls + 1;
while exitflag == 1 && fval > tol % Feasible and some slacks are nonzero
    c = 0;
    for i = 1:mi        
        j = ineq_slack_offset+i;
        if x(j) > tol
            ub_r(j) = 0.0;
            ineq_iis(i) = true;
            c = c+1;
        end
    end
    for i = 1:me
        j = eq_pos_slack_offset+i;
        if x(j) > tol            
            ub_r(j) = 0.0;
            eq_iis(i) = true;
            c = c+1;
        end
    end
    for i = 1:me
        j = eq_neg_slack_offset+i;
        if x(j) > tol            
            ub_r(j) = 0.0;
            eq_iis(i) = true;
            c = c+1;
        end
    end
    [x,fval,exitflag] = linprog(f,Aineq_r,bineq,Aeq_r,beq,lb_r,ub_r,[],opts);
    if fval > 0
        fval0 = fval;
    end
    ncalls = ncalls + 1;
end 
end

Этот код создает generatecover функция помощника. Код использует тот же метод индексирования для отслеживания ограничений, что и удаление IIS в коде цикла.

function [coversetineq,coverseteq,nlp] = generatecover(Aineq,bineq,Aeq,beq,lb,ub)
% Returns the cover set of linear inequalities, the cover set of linear
% equalities, and the total number of calls to linprog.
% Adapted from Chinneck [1] Algorithm 7.3. Step numbers are from this book.
coversetineq = [];
coverseteq = [];
activeA = true(size(bineq));
activeAeq = true(size(beq));
% Step 1 of Algorithm 7.3
[ineq_iis,eq_iis,ncalls] = elasticfilter(Aineq,bineq,Aeq,beq,lb,ub);
nlp = ncalls;
ninf = sum(ineq_iis(:)) + sum(eq_iis(:));
if ninf == 1
    coversetineq = ineq_iis;
    coverseteq = eq_iis;
    return
end
holdsetineq = find(ineq_iis);
holdseteq = find(eq_iis);
candidateineq = holdsetineq;
candidateeq = holdseteq;
% Step 2 of Algorithm 7.3
while sum(candidateineq(:)) + sum(candidateeq(:)) > 0
    minsinf = inf;
    ineqflag = 0;
    for i = 1:length(candidateineq(:))
        activeA(candidateineq(i)) = false;
        idx2 = find(activeA);
        idx2eq = find(activeAeq);
        [ineq_iis,eq_iis,ncalls,fval] = elasticfilter(Aineq(activeA,:),bineq(activeA),Aeq(activeAeq),beq(activeAeq),lb,ub);
        nlp = nlp + ncalls;
        ineq_iis = idx2(find(ineq_iis));
        eq_iis = idx2eq(find(eq_iis));
        if fval == 0
            coversetineq = [coversetineq;candidateineq(i)];
            return
        end
        if fval < minsinf
            ineqflag = 1;
            winner = candidateineq(i);
            minsinf = fval;
            holdsetineq = ineq_iis;
            if numel(ineq_iis(:)) + numel(eq_iis(:)) == 1
                nextwinner = ineq_iis;
                nextwinner2 = eq_iis;
                nextwinner = [nextwinnner,nextwinner2];
            else
                nextwinner = [];
            end
        end
        activeA(candidateineq(i)) = true;
    end
    for i = 1:length(candidateeq(:))
        activeAeq(candidateeq(i)) = false;
        idx2 = find(activeA);
        idx2eq = find(activeAeq);
        [ineq_iis,eq_iis,ncalls,fval] = elasticfilter(Aineq(activeA),bineq(activeA),Aeq(activeAeq,:),beq(activeAeq),lb,ub);
        nlp = nlp + ncalls;
        ineq_iis = idx2(find(ineq_iis));
        eq_iis = idx2eq(find(eq_iis));
        if fval == 0
            coverseteq = [coverseteq;candidateeq(i)];
            return
        end
        if fval < minsinf
            ineqflag = -1;
            winner = candidateeq(i);
            minsinf = fval;
            holdseteq = eq_iis;
            if numel(ineq_iis(:)) + numel(eq_iis(:)) == 1
                nextwinner = ineq_iis;
                nextwinner2 = eq_iis;
                nextwinner = [nextwinnner,nextwinner2];
            else
                nextwinner = [];
            end
        end
        activeAeq(candidateeq(i)) = true;
    end
% Step 3 of Algorithm 7.3
    if ineqflag == 1
        coversetineq = [coversetineq;winner];
        activeA(winner) = false;
        if nextwinner
            coversetineq = [coversetineq;nextwinner];
            return
        end
    end
    if ineqflag == -1
        coverseteq = [coverseteq;winner];
        activeAeq(winner) = false;
        if nextwinner
            coverseteq = [coverseteq;nextwinner];
            return
        end
    end
    candidateineq = holdsetineq;
    candidateeq = holdseteq;
end
end

См. также

linprog

Связанные темы

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

Документация