В этом примере показано, как исследовать линейные ограничения, которые вызывают недопустимость задачи. Для получения дополнительной информации об этих методах см. 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)
, где ограничения равенства имеют как аддитивные, так и вычитающие положительные упругие значения.
Затем решите связанную задачу линейного программирования (LP)
удовлетворяющее перечисленным ограничениям и с .
Если у связанной НД есть решение, удалите все ограничения, которые строго положительно связаны , и записать эти ограничения в список индексов (потенциальные представители IIS). Вернитесь к предыдущему шагу, чтобы решить новую, уменьшенную связанную НД.
Если связанная НД не имеет решения (недопустимо) или не связана строго положительная связь , выход из фильтра.
Эластичный фильтр может выйти во многих меньших итерациях, чем фильтр удаления, потому что он может принести в 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 за раз из задачи, пока проблема не станет допустимой. Код использует метод индексации, чтобы отслеживать ограничения с точки зрения их позиций в исходной задаче, прежде чем алгоритм удаляет какие-либо ограничения.
Код отслеживает исходные переменные в задаче с помощью логического вектора 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 непосредственно. Как объясняет 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