Исследуйте линейный Infeasibilities

Скрипт Open Live Script

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

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

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

Неосуществимый пример

Создайте случайный матричный 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] = 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 линейных ограничений, таким образом, функция под названием linprog 150 раз.

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

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

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

$\begin{array}{l} A_{i n e q} x \leq b_{i n e q} + e \\ A_{e q} x = b_{e q} + e_{1} - e_{2} \end{array}$

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

$\min_{x, e} \sum e_{i}$

подвергните перечисленным ограничениям и с $e_{i} \geq 0$ .

Если связанный LP имеет решение, удалите все ограничения, которые имеют сопоставленное строго положительное $e_{i}$ , и запишите те ограничения в списке индексов (потенциальные члены IIS). Возвратитесь к предыдущему шагу, чтобы решить новое, уменьшал сопоставленный LP.
Если связанный LP не имеет никакого решения (неосуществимо), или имеет не строго положительный сопоставленный $e_{i}$ , выйдите из фильтра.

Эластичный фильтр может выйти во многих из меньшего количества итераций, чем фильтр удаления, потому что это может принести много индексов целиком в 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)

Удалите 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, чтобы найти набор покрытия ограничений, который, когда удалено, дает выполнимый набор. Алгоритм реализован в generatecover функция помощника в конце этого примера.

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

coversetineq = 5×1

   114
    97
    60
    82
     2

coverseteq =

     []

nlp = 40

Удалите эти ограничения и решите LP.

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] 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 функция помощника. Код использует тот же метод индексации в отслеживании ограничения как Удалить 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

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