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

Открыть Live Script

В этом примере показано, как исследовать линейные ограничения, которые вызывают недопустимость задачи. Для получения дополнительной информации об этих методах см. 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 вспомогательная функция в конце этого примера.

Примечание: Если вы используете файл live скрипта для этого примера, 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

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

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

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

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

Aineqx≤bineq+eAeqx=beq+e1-e2

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

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, является IIS. Этот результат согласуется с ответом фильтра удаления. Различие между подходами заключается в том, что комбинированный подход эластичного и удаляющего фильтра использует намного меньше linprog вызовы. Отображение общего количества linprog вызовы, используемые эластичным фильтром, за которыми следует фильтр удаления.

disp(ncall + ncall2)
     7

Удаление 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 непосредственно. Как объясняет Chinneck [1], нахождение MaxFS является NP-полной задачей, что означает, что задача не обязательно имеет эффективные алгоритмы для нахождения MaxFS. Однако Чиннек предлагает некоторые алгоритмы, которые могут работать эффективно.

Используйте Алгоритм 7.3 Chinneck, чтобы найти набор ограничений, который при удалении дает допустимый набор. Алгоритм реализован в 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 Chinneck начинается с результата упругого фильтра и затем сохраняет только необходимые ограничения.

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

Кроме того, заметьте, что неравенства, удаленные в цикле, не совсем совпадают с неравенствами, удаленными Алгоритмом 7.3. Цикл удаляет неравенства 114, 97, 82, 60 и 64, в то время как Алгоритм 7.3 удаляет неравенства 114, 97, 82, 60 и 2. Проверяйте, что неравенства 82 и 64 образуют IIS (как указано в Remove 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] Chinneck, J. W. Выполнимость и недопустимость в оптимизации: алгоритмы и вычислительные методы. Спрингер, 2008.

[2] Chinneck, J. W. «Осуществимость и недопустимость в оптимизации». Руководство по 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 вспомогательная функция. Код использует тот же метод индексации для отслеживания ограничений, что и Remove 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

См. также

Похожие темы

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

Поддержка

Для просмотра документации необходимо авторизоваться на сайте