혼합 정수 선형 계획법(MILP) 알고리즘
혼합 정수 선형 계획법 정의
혼합 정수 선형 계획법(MILP)은 다음 조건을 갖는 문제입니다.
선형 목적 함수 fTx. 여기서 f는 상수로 구성된 열 벡터이고 x는 미지수로 구성된 열 벡터입니다.
범위와 선형 제약 조건. 단, 비선형 제약 조건은 포함되지 않음(정의는 제약 조건 작성하기 참조)
x의 일부 성분이 정수 값을 가져야 한다는 제한
수학적 측면에서 벡터 f, lb, ub와 행렬 A, Aeq, 이에 대응하는 벡터 b, beq, 그리고 인덱스 세트 intcon이 주어질 때, 다음 문제의 해가 되는 벡터 x를 구합니다.
HiGHS MILP 알고리즘
HiGHS 개요
intlinprog 알고리즘은 HiGHS 오픈소스 소프트웨어를 기반으로 합니다. intlinprog는 MATLAB® 형식의 입력값과 옵션을 동등한 HiGHS 인수로 변환하고, 반환된 해도 표준 MATLAB 형식으로 변환합니다.
알고리즘 개요
알고리즘은 다음 단계를 수행합니다.
분기/한정 노드 가져오기(없는 경우 중지)
다음 중지 조건까지 반복:
다음 중지 조건까지 탐색:
도메인 전파
노드를 가지치기하고, 범위를 업데이트하고, 실현불가능성을 감지하거나
ObjectiveCutOff옵션 값에 도달한 경우 종료재시작 조건을 충족하는 경우 2단계로 돌아감
다음 노드 설치:
노드 선택 및 평가
실행에서 이 노드를 가지치기하는 경우, 5단계로 돌아감
노드에 대한 절단 생성
도메인이 실현 가능하지 않은 경우 노드를 자르고 열린 노드들을 노드 큐로 가져옴
기저 업데이트
4단계로 이동
풀이 전처리
일반적으로, 문제에 포함되는 변수의 개수(x의 성분 개수)를 줄이고 선형 제약 조건의 개수를 줄일 수 있습니다. 이렇게 줄이는 작업을 수행하느라 솔버에서 시간이 더 걸릴 수 있지만, 일반적으로 해를 구하는 데 소요되는 전체 시간이 감소되므로 더 큰 문제를 풀 수 있게 됩니다. 알고리즘은 해를 수치적으로 더 안정하게 만들 수 있습니다. 또한, 이러한 알고리즘은 경우에 따라 실현불가능 문제를 감지하기도 합니다.
풀이 전처리 단계는 중복된 변수와 제약 조건을 제거하고, 모델의 스케일링과 제약 조건 행렬의 희소성을 향상시키고, 변수의 범위를 강화하고, 모델의 원문제 실현불가능성과 쌍대 문제 실현불가능성을 감지하는 것을 목표로 합니다. 배경 정보는 Andersen과 Andersen의 문헌[3] 및 Mészáros와 Suhl의 문헌[10]을 참조하십시오.
혼합 정수 계획의 전처리 과정에서 intlinprog는 정수 제한과 함께 선형 부등식 A*x ≤ b를 분석하여 다음을 결정합니다.
문제가 실현 불가능한지 여부.
일부 범위를 좁힐 수 있는지 여부.
일부 부등식이 중복되어 있어, 무시하거나 제거할 수 있는지 여부.
일부 부등식을 강화할 수 있는지 여부.
일부 정수 변수를 고정할 수 있는지 여부.
정수 전처리에 대한 배경 정보는 Achterberg 등이 작성한 문헌[1]을 참조하십시오.
루트 노드 평가하기
루트 노드를 평가하기 위해 알고리즘은 다음 단계를 수행합니다.
대칭을 감지하고 문제를 단순화합니다.
루트 LP를 평가합니다.
LP 절단을 생성하고 추가합니다(절단 생성 참조).
무작위 반올림을 적용합니다.
LP 절단을 더 생성하고 추가합니다.
루프에서 절단 생성 및 발견법을 수행합니다.
재시작 조건을 확인하고 필요한 경우 2단계에서 재시작합니다.
루트 노드 평가를 완료하면 알고리즘은 분기한정으로 넘어갑니다.
절단 생성
절단은 intlinprog가 문제에 추가하는 부가적인 선형 부등식 제약 조건입니다. 이 부등식은 해가 정수에 더 가깝도록 LP 완화의 실현 가능한 영역을 제한하려고 시도합니다. 절단 생성 알고리즘(평면 절단 방법이라고도 함)에 대한 배경 정보는 Cornuéjols의 문헌[6] 및 Atamtürk, Nemhauser, Savelsbergh의 문헌[4]을 참조하십시오.
플런지/급강하
정수 실현가능점을 구하기 위해 intlinprog는 분기한정 스텝과 유사한 발견법을 사용하되, 다른 분기를 만들지 않고 트리에서 바로 아래에 있는 분기를 따릅니다. 이 단일 분기는 트리 조각 아래로 신속한 “급강하”를 초래하기 때문에 “급강하(Diving)”라는 이름이 붙었습니다.
급강하 발견법은 일반적으로 정수 값을 가져야 하는 하나의 변수를 선택합니다. 이에 대한 현재 해는 소수입니다. 그런 다음 이 발견법에서는 변수가 정수 값을 갖도록 강제하는 범위를 추가하고 연관된 완화된 LP를 다시 풉니다. 급강하 발견법 간의 주요 차이점은 범위에 대한 변수를 선택하는 방법에 있습니다. Berthold의 문헌[5], 섹션 3.1을 참조하십시오.
무작위 반올림, RINS, RENS
새로운 정수 실현가능점을 구하기 위해 intlinprog는 현재 최적의 실현 가능한 정수 해에 해당하는 점(사용 가능한 경우)의 근방을 탐색하여 새롭고 더 나은 해를 구합니다. Danna, Rothberg, Le Pape의 문헌[8]을 참조하십시오. 유사하게 새로운 정수 실현가능점을 구하기 위해 intlinprog는 한 노드에서 완화된 문제의 LP 해를 얻고, 실현가능성을 유지하기 위한 시도로서 정수 성분을 반올림합니다. 무작위 반올림 스텝을 통해 intlinprog는 경우에 따라 새로운 실현가능점을 구할 수 있습니다. RENS(Relaxation Enforced Neighborhood Search)는 정수 실현가능점을 찾는 또 다른 탐색 기법입니다. Berthold의 문헌[5]을 참조하십시오.
분기한정
분기한정법은 MILP의 해에 수렴하려고 시도하는 일련의 하위 문제를 생성합니다. 하위 문제는 해 fTx에 대한 일련의 상한 및 하한을 제공합니다. 이러한 범위를 원문제 범위와 쌍대 문제 범위라고 부릅니다. 최소화 문제의 경우 첫 번째 상한(원문제)은 임의의 실현 가능한 해이고, 첫 번째 하한(쌍대 문제)은 완화된 문제의 해입니다. 최대화 문제의 경우 원문제 범위는 하한이고 쌍대 문제 범위는 상한입니다. 완화된 선형 계획법 문제의 해는 MILP의 해보다 낮은 목적 함수 값을 가집니다. 또한, 임의의 실현가능점 xfeas는 다음을 충족합니다.
| fTxfeas ≥ fTx, | (1) |
그 이유는 fTx가 모든 실현가능점 중 최솟값이기 때문입니다.
이 컨텍스트에서 노드는 원래 문제와 동일한 목적 함수, 범위, 선형 제약 조건을 갖지만 정수 제약 조건은 갖지 않으며 선형 제약 조건 또는 범위에 대한 특정한 변경 사항이 적용된 LP입니다. 루트 노드는 정수 제약 조건이 없고 선형 제약 조건 또는 범위에 대한 변경 사항이 적용되지 않은 원래 문제입니다. 즉, 루트 노드는 완화된 초기 LP입니다.
시작 범위에서 분기한정법은 루트 노드에서 분기를 생성함으로써 새 하위 문제를 생성합니다. 분기 생성 단계는 여러 규칙 중 하나에 따라 발견법 방식으로 수행됩니다. 각 규칙은 한 변수가 정수 J보다 작거나 같도록 제한하거나 J+1보다 크거나 같도록 제한하여 문제를 분할하는 아이디어를 기반으로 합니다. intcon에 지정된 정수에 대응하는 x LP의 요소가 정수가 아닌 경우 이 두 하위 문제가 생성됩니다. 여기서, xLP는 완화된 문제의 해입니다. 변수를 내림(버림)한 값을 J로 취하고 올림(반올림)한 값을 J+1로 취합니다. 결과로 생성되는 두 문제는 fTxLP보다 크거나 같은 해를 가집니다. 그 이유는 제한 사항이 더 많기 때문입니다. 따라서 최소화 문제의 경우 이 절차는 잠재적으로 하한을 높일 수 있습니다.
알고리즘이 분기를 생성하고 나면 탐색할 노드 두 개가 새로 생성됩니다. intlinprog는 목적 함수 값과 현재 전역 범위를 비교하여 일부 하위 문제의 분석을 건너뜁니다.
분기한정 절차는 다음 중지 기준 중 하나가 충족될 때까지 분석할 하위 문제를 체계적으로 생성하고 목적 함수의 상한 또는 하한을 향상시키지 않는 하위 문제를 삭제하는 작업을 계속 수행합니다.
문제가 실현 불가능한 것으로 증명됩니다.
목적 함수 값이
ObjectiveCutOff제한에 도달합니다.알고리즘이
MaxTime옵션을 초과합니다.목적 함수의 하한과 상한 사이의 차이가
AbsoluteGapTolerance허용오차 또는RelativeGapTolerance허용오차보다 작습니다.탐색된 노드 개수가
MaxNodes옵션을 초과합니다.
분기한정 절차에 대한 배경 정보는 Nemhauser와 Wolsey의 문헌[11] 및 Wolsey의 문헌[13]을 참조하십시오.
반복 과정 표시
Display 옵션을 디폴트 "iter"으로 설정하여 반복 과정 표시를 선택할 경우 솔버는 일부 단계를 표시합니다. HiGHS의 반복 과정 표시는 기타 솔버의 반복 과정 표시보다 더 광범위하고 복잡합니다. 또한 HiGHS 알고리즘은 분기한정 탐색을 다시 시작하여 반복 과정 표시도 다시 시작되게 할 수 있습니다.
반복 과정 표시를 선택하는 방법:
options = optimoptions("intlinprog",Display="iter"); [x,fval,exitflag,output] = intlinprog(f,intcon,A,b,Aeq,beq,... lb,ub,options)
프리앰블. 반복 과정 표시는 "풀이 전처리"의 결과를 표시하는 것에서 시작합니다. 풀이 전처리 알고리즘은 선형 제약 조건 행렬에서 중복되는 행과 열을 식별하고 제거하며 문제와 관련된 단순화를 수행하여 원래 문제의 복잡도를 줄입니다. 예를 들면 다음과 같습니다.
Presolving model 18018 rows, 26027 cols, 248579 nonzeros 15092 rows, 24343 cols, 217277 nonzeros Objective function is integral with scale 1
루트 노드 평가. 루트 노드는 정수 제약 조건을 고려하지 않은, 문제의 선형 계획법 해입니다. 최소화 문제의 경우 루트 노드의 목적 함수 값은 정수 제약 조건을 포함하는 문제에 대한 해의 목적 함수 값의 하한입니다. 최소화 문제의 경우 상한(있는 경우)은 모든 제약 조건에 대한 실현가능점에서 가져옵니다. 발견된 실현가능점이 아직 없는 경우 상한은 Inf입니다.
반복 과정 표시는 풀이 전처리 후 문제의 크기를 보여줍니다.
Solving MIP model with: 15092 rows 24343 cols (24343 binary, 0 integer, 0 implied int., 0 continuous) 217277 nonzeros
binary는 이진 변수의 개수입니다.integer는 정수 변수의 개수입니다.implied int는 정수로 유추되는 변수 개수입니다. 예를 들어,x(1)이 정수이고x(1) + x(2) = 5인 경우,x(2)는 정수로 유추됩니다.continuous는 연속 변수의 개수입니다.
변수의 총 개수는 모델에 있는 열의 개수입니다.
동적 제약 조건 생성. 소프트웨어는 반복 과정 표시에 다음과 같은 세 개의 헤더가 있는 "동적 제약 조건"의 생성으로 시작합니다.
Cuts — 활성 절단 개수
inLp — LP 행렬에서 비모델 행의 개수
Confl. — 충돌 개수
소프트웨어는 동적 제약 조건 생성을 다시 시작하여 제약 조건을 더욱 확장할 수 있으며, 이는 생성 프로세스를 새 상태에서 시작하여 더 많은 제약 조건을 생성할 수 있습니다.
분기한정 과정을 시작하기 전의 마지막 스텝에서 소프트웨어는 대칭 감지 결과와 발견된 생성기 및 오비토프(orbitope)의 개수를 보고합니다. 예를 들어, 다음은 프리앰블 및 동적 제약 조건 생성 단계에 나타나는 반복 과정 표시의 일부입니다.
Nodes | B&B Tree | Objective Bounds | Dynamic Constraints | Work
Proc. InQueue | Leaves Expl. | BestBound BestSol Gap | Cuts InLp Confl. | LpIters Time
0 0 0 0.00% 0 inf inf 0 0 0 0 1.5s
0 0 0 0.00% 0 inf inf 0 0 5 4934 3.6s
…
0 0 0 0.00% 0 inf inf 4630 553 289 140296 159.6s
0.0% inactive integer columns, restarting
Model after restart has 15090 rows, 24338 cols (24338 bin., 0 int., 0 impl., 0 cont.), and 217075 nonzeros
0 0 0 0.00% 0 inf inf 550 0 0 140296 160.5s
…
0 0 0 0.00% 0 inf inf 5602 524 260 277323 318.0s
6.2% inactive integer columns, restarting
Model after restart has 14185 rows, 22816 cols (22816 bin., 0 int., 0 impl., 0 cont.), and 200423 nonzeros
0 0 0 0.00% 0 inf inf 524 0 0 277323 318.6s
…
0 0 0 0.00% 0 inf inf 4683 535 525 408393 505.8s
Symmetry detection completed in 3.1s
Found 215 generators and 12 full orbitope(s) acting on 730 columns대칭 감지와 오비토프에 대한 정보는 Hojny, Pfetsch, Schmitt의 문헌[7] 및 Pfetsch와 Rehn의 문헌[12]을 참조하십시오. 소프트웨어는 다음에 설명된 것처럼 분기한정 스텝을 진행하면서 동적 제약 조건을 계속 추가합니다.
분기한정. 분기한정법은 MILP의 해에 수렴하려고 시도하는 일련의 하위 문제를 생성합니다. 하위 문제는 해 fTx에 대한 일련의 상한 및 하한을 제공합니다. 최소화 문제의 경우 첫 번째 상한은 임의의 실현 가능한 해이고, 첫 번째 하한은 완화된 문제의 해입니다.
분기한정 절차 도중에 intlinprog는 다음 반복 과정 표시를 생성합니다.
Nodes | B&B Tree | Objective Bounds | Dynamic Constraints | Work
Proc. InQueue | Leaves Expl. | BestBound BestSol Gap | Cuts InLp Confl. | LpIters Time
72 0 2 1.56% 0 inf inf 4695 535 738 667009 686.9s
…
T 271 107 51 1.56% 0 449 100.00% 6051 271 1767 792786 776.8s
T 279 107 53 1.56% 0 439 100.00% 6053 271 1794 793241 777.5s
…
L 1223 538 295 1.98% 0 434 100.00% 6984 243 7628 1689k 1580.6s
…
1321 539 333 1.98% 0 434 100.00% 7029 243 8628 1898k 1650.6s
Restarting search from the root node
Model after restart has 13947 rows, 22426 cols (22426 bin., 0 int., 0 impl., 0 cont.), and 194633 nonzeros
1323 0 0 0.00% 0 434 100.00% 243 0 0 1902k 1653.5s
1323 0 0 0.00% 0 434 100.00% 243 76 10 1905k 1655.7s
…
1694 173 52 0.00% 0 434 100.00% 9411 318 2220 3205k 2584.3s
B 1710 167 55 0.00% 0 433 100.00% 9415 318 2263 3207k 2586.2s
1726 224 56 0.00% 0 433 100.00% 9999 378 2307 3237k 2608.7s가장 왼쪽 열에는 표시의 해당 행의 새로운 실현가능점을 찾은 방법을 나타내는 다음과 같은 코드가 표시됩니다.
L— 원문제 발견법 도중 하위 MIP 문제를 푸는 동안T— 노드를 평가하는 동안, 트리 탐색 도중에B— 분기 생성 도중H— 발견법을 사용하여P— MIP를 풀기 전 시작 도중C— 중심 반올림을 사용하여R— 무작위 반올림을 사용하여S— LP를 푸는 동안F— 실현가능성 펌프 사용U— 비유계 LP에서
마침. 알고리즘이 중지된 이유와 결과 요약이 반복 과정 표시의 끝부분에 나타납니다.
Solving report
Status Time limit reached
Primal bound 14
Dual bound 0
Gap 100% (tolerance: 0.01%)
Solution status feasible
14 (objective)
0 (bound viol.)
1.7763568394e-15 (int. viol.)
0 (row viol.)
Timing 7200.02 (total)
3.08 (presolve)
0.00 (postsolve)
Nodes 5830
LP iterations 9963775 (total)
2657448 (strong br.)
324908 (separation)
2140011 (heuristics)
Solver stopped prematurely. Integer feasible point found.
Intlinprog stopped because it exceeded the time limit, options.MaxTime = 7200 . The intcon variables are integer within tolerance,
options.ConstraintTolerance = 1e-06.Status— 반복이 중지된 이유Primal bound— 최소화 문제에 대한 목적 함수 값의 상한, 최대화 문제에 대한 목적 함수 값의 하한Dual bound— 최대화 문제에 대한 목적 함수 값의 하한, 최소화 문제에 대한 목적 함수 값의 상한Gap— 원문제 범위와 쌍대 문제 범위 사이의 상대적 간격Solution status해의 상태
목적 함수 값으로,
(objective)로 레이블이 지정됨변수 범위에 대한 해의 최대 위반으로,
(bound viol.)로 레이블이 지정됨정수 값에서 정수형 변수의 최대 위반으로,
(int. viol.)로 레이블이 지정됨선형 제약 조건에 대한 해의 최대 위반으로,
(row viol.)로 레이블이 지정됨
Timing— 다양한 해 단계의 시간 설정(단위: 초)Nodes— 탐색된 노드 개수LP iterations— 단계별 선형 계획 반복 횟수와 총 선형 계획 반복 횟수
참고 문헌
[1] Achterberg, T.,Bixby, R., Gu, Z., Rothberg, E., and Weninger, D. Presolve Reductions in Mixed Integer Programming. INFORMS J. Computing 32(2), 2019, pp. 473–506. Available at https://doi.org/10.1287/ijoc.2018.0857.
[2] Achterberg, T., T. Koch and A. Martin. Branching rules revisited. Operations Research Letters 33, 2005, pp. 42–54. Available at https://opus4.kobv.de/opus4-zib/files/788/ZR-04-13.pdf.
[3] Andersen, E. D., and Andersen, K. D. Presolving in linear programming. Mathematical Programming 71, pp. 221–245, 1995.
[4] Atamtürk, A., G. L. Nemhauser, M. W. P. Savelsbergh. Conflict graphs in solving integer programming problems. European Journal of Operational Research 121, 2000, pp. 40–55.
[5] Berthold, T. Primal Heuristics for Mixed Integer Programs. Technischen Universität Berlin, September 2006. Available at https://opus4.kobv.de/opus4-zib/files/1029/Berthold_Primal_Heuristics_For_Mixed_Integer_Programs.pdf.
[6] Cornuéjols, G. Valid inequalities for mixed integer linear programs. Mathematical Programming B, Vol. 112, pp. 3–44, 2008.
[7] Hojny, C., Pfetsch, M. E., Schmitt, A. Extended formulations for column-constrained orbitopes. TU Darmstadt, Department of Mathematics, Research Group Optimization, Dolivostr. 15, 64293 Darmstadt, June 2017. Available at https://optimization-online.org/2017/06/6092/.
[8] Danna, E., Rothberg, E., Le Pape, C. Exploring relaxation induced neighborhoods to improve MIP solutions. Mathematical Programming, Vol. 102, issue 1, pp. 71–90, 2005.
[9] Hendel, G. New Rounding and Propagation Heuristics for Mixed Integer Programming. Bachelor's thesis at Technische Universität Berlin, 2011. PDF available at https://opus4.kobv.de/opus4-zib/files/1332/bachelor_thesis_main.pdf.
[10] Mészáros C., and Suhl, U. H. Advanced preprocessing techniques for linear and quadratic programming. OR Spectrum, 25(4), pp. 575–595, 2003.
[11] Nemhauser, G. L. and Wolsey, L. A. Integer and Combinatorial Optimization. Wiley-Interscience, New York, 1999.
[12] Pfetsch, M. E. and Rehn, T. A computational comparison of symmetry handling methods for mixed integer programs. Mathematical Programming Computation (2019) 11:37–93. Available at https://optimization-online.org/wp-content/uploads/2015/11/5209.pdf.
[13] Wolsey, L. A. Integer Programming. Wiley-Interscience, New York, 1998.
