含分式的非线性规划求解: 基于Outer Approximation的Branch-and-cut 算法及其Java实现_outer approximation algorithm

问题介绍

在混合整数规划求解中，如果目标函数和约束都为线性的情况下，求解器可以直接求解得到精确解。虽然目前求解器也可以完成一些特殊非线性目标函数和约束的求解，如MIQP, MIQCP, MIQCQP, MISOCP等（详情参见运小筹第60期推文），但针对一般的非线性问题，如目标函数为以下的分式形式，就无法利用求解器直接求解。
$$min\sum _{i\in I}\frac{U_1}{U_2+{\sum }_{j\in J}{x}_j}\text{\hspace{1em}(1)}$$
针对目标函数为凸函数的情况，我们采取Outer-Approximation(下称OA)的方法进行求解，OA由Duran和Grossmann在1986年提出${}^{[1]}$，后面被不断发展改进。OA的基本思想就是利用松弛问题的解生成切线，再加入松弛主问题中，直至找到最优解，利用多条切线去刻画原来的凸函数曲线。近年来的研究中，学者们通常将OA与Branch-and-cut结合使用，利用OA生成cut加入算法中。

Outer Approximation介绍

下面将借助一个具体模型介绍OA的求解过程。
$$min\sum _{i\in I}\frac{U_1}{U_2+{\sum }_{j\in J}{\omega }_{ij}{x}_j}\text{\hspace{1em}(2)}$$
$$\sum _{j\in J}{x}_j\le r\text{\hspace{1em}(3)}$$
x j ∈ { 0 , 1 } , ∀ j ∈ J (4) x_j \in \{0,1\}, \forall j\in J\tag{4} xj​∈{0,1},∀j∈J(4)
针对原问题，我们进行一些变换，引入连续变量$\theta$，将复杂的目标函数放在约束条件中
$$min\theta \text{\hspace{1em}(5)}$$
$$\theta \ge \sum _{i\in I}\frac{U_1}{U_2+{\sum }_{j\in J}{\omega }_{ij}{x}_j}\text{\hspace{1em}(6)}$$
$$\sum _{j\in J}{x}_j\le r\text{\hspace{1em}(7)}$$
x j ∈ { 0 , 1 } , ∀ j ∈ J (8) x_j \in \{0,1\}, \forall j\in J \tag{8} xj​∈{0,1},∀j∈J(8)
对约束(6)，它在$x$上是凸函数，所以我们在$x^{\ast }$处进行一阶展开得到
$$\theta \ge F(x^{\ast })+\sum _{j\in J}{F}_j^{{}^{\prime }}(x^{\ast })(x_j-x_j^{\ast })\text{\hspace{1em}(9)}$$
其中$F(x)={\sum }_{i\in I}\frac{U_1}{U_2+{\sum }_{j\in J}{\omega }_{ij}{x}_j}$
$$F_j^{{}^{\prime }}(x)=\frac{dF(x)}{d{x}_j}=-\sum _{i\in I}\frac{U_1{\omega }_{ij}}{(U_2+{\sum }_{k\in J}{w}_{ik}{x}_k{)}^2}\text{\hspace{1em}(10)}$$
所以约束(6)变换为
$$\theta \ge \sum _{i\in I}\frac{U_1}{U_2+{\sum }_{j\in J}{\omega }_{ij}{x}_j^{\ast }}-\sum _{j\in J}\sum _{i\in I}\frac{U_1{\omega }_{ij}}{(U_2+{\sum }_{k\in J}{w}_{ik}{x}_k^{\ast }{)}^2}(x_j-x_j^{\ast })\text{\hspace{1em}(11)}$$
式(11)就被称为OA cut

java调用cplex代码

Branch-and-cut实现: 为了求解上述模型，我们借助cplex求解器中的callback函数，当当前LP松弛的解决方案$x^{\ast }$被证明是整数时，我们检查是否违反了约束(11)，如果违背，则将它们添加到当前LP中，OA cut是全局有效的，通过在MILP求解器的lazy-cut callback实现的。下面是Java调用cplex实现OA cut的Branch-and-cut算法代码。
首先是Cut类，用于存储OA生成的cuts

import java.util.ArrayList;

import ilog.concert.IloNumExpr;

public class Cuts {
	ArrayList<IloNumExpr> lhs;
	ArrayList<Double> rhs;
	/**
	 * @return the lhs
	 */
	public ArrayList<IloNumExpr> getLhs() {
		return lhs;
	}
	/**
	 * @param lhs the lhs to set
	 */
	public void setLhs(ArrayList<IloNumExpr> lhs) {
		this.lhs = lhs;
	}
	/**
	 * @return the rhs
	 */
	public ArrayList<Double> getRhs() {
		return rhs;
	}
	/**
	 * @param rhs the rhs to set
	 */
	public void setRhs(ArrayList<Double> rhs) {
		this.rhs = rhs;
	}

}

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35

下面是主要的代码

import java.io.IOException;
import java.util.ArrayList;
import java.util.Arrays;
import Cuts;
import ilog.concert.*;
import ilog.cplex.IloCplex;
import ilog.cplex.IloCplex.LazyConstraintCallback;

public class OA {
	
	//calculate the gradient
	public static double CalculateGradient(double[] x,double[][] w,int j,int num,double U1,double U2) throws IloException {
		double gradient_j=0; 
		for(int i=0;i<num;i++) {
			double b=0;
			for(int k=0;k<x.length;k++) {
				b+=w[i][k]*x[k];
			}
			
			gradient_j+=-U1*w[i][j]/(b+U2)/(b+U2);
			
		}
		return gradient_j;
	}
	
	
	//calculate the value of the objective function
	public static double calF(double[] x,double[][] w,int num,double U1,double U2) {
		double F=0;
		for(int i=0;i<num;i++) {
			double b=0;
			for(int k=0;k<x.length;k++) {
				b+=w[i][k]*x[k];
			}
			
			F+= U1/(U2+b);
			
		}
		return F;
	}
	//OA cut
	public static Cuts makecuts(IloNumVar[] x, double[] xSol,double[][] w,int num,double U1,double U2, IloModeler ilcplex,IloNumVar theta,double theta0) throws IloException {
		Cuts cuts=new Cuts();
		ArrayList<IloNumExpr> cutLhs = new ArrayList<IloNumExpr>();
		ArrayList<Double> cutRhs = new ArrayList<Double>();
		
		IloLinearNumExpr lhsExpr=ilcplex.linearNumExpr();
		double rhsExpr=0;
		double lhs = 0, rhs = 0;
		
		
		lhsExpr.addTerm(1, theta); 
		
		rhsExpr+=calF(xSol,w,num,U1,U2);
		
		
		for(int j=0;j<xSol.length;j++) {
			lhsExpr.addTerm(-CalculateGradient(xSol,w,j,num,U1,U2), x[j]);
			rhsExpr-=CalculateGradient(xSol,w,j,num,U1,U2)*xSol[j];
		}
		lhs=theta0;
		rhs=calF(xSol,w,num,U1,U2);

		if(lhs<rhs) {
			cutLhs.add(lhsExpr);
			cutRhs.add(rhsExpr);
		}
		
		cuts.setLhs(cutLhs);
		cuts.setRhs(cutRhs);
		
		return cuts;
	}

	public static void buildModel(double[][] w,int num1,int num2,double U1,double U2) throws IloException {
		IloCplex ilcplex =new IloCplex();
	    IloIntVar[] x=new IloIntVar[num1];
	    int r=3;

	    for(int i=0;i<num1;i++) {
	    	x[i]=ilcplex.boolVar();
	    }
	    
	    
	    
	  //objective
	    IloNumVar theta=ilcplex.numVar(0, Double.MAX_VALUE);
	    
	    ilcplex.addMinimize(theta);
	  //constraints
	    IloLinearNumExpr left=ilcplex.linearNumExpr();
	    for(int i=0;i<num1;i++) {
	    	left.addTerm(1, x[i]);
	    }
	    ilcplex.addEq(left, r);
	    
 		ilcplex.use(new LazyCallback(x,ilcplex,theta,w,num2,U1,U2));
 		
 		
	    
	    if(ilcplex.solve()) {
	    	double[] xVal = ilcplex.getValues(x);
	    	
	    	System.out.println("The objective value is:" + ilcplex.getObjValue());
			
			System.out.println(Arrays.toString(xVal));
			
			
			
	    }
	}
	//LazyCallback
	public static class LazyCallback extends LazyConstraintCallback{
		Cuts cut;
		ArrayList<IloNumExpr> cutLhs;
		ArrayList<Double> cutRhs;
		IloIntVar[] x;
		IloCplex ilcplex;
		IloNumVar theta;
		double[][] w;
		int num;
		double U1,U2;
		
		
		
		LazyCallback(IloIntVar[] x0,IloCplex ilcplex0,IloNumVar theta0,double[][]w0,int num0,double U10,double U20){
			x=x0;
			ilcplex=ilcplex0;
			theta=theta0;
			w=w0;
			num=num0;
			U1=U10;
			U2=U20;
			
		}
		
		public void main() throws IloException{
			double[] xSol =getValues(x);
			double theta0=getValue(theta);
			cut = makecuts(x, xSol,w,num,U1,U2,ilcplex,theta,theta0);
		
			cutLhs = cut.getLhs();
			cutRhs= cut.getRhs();
			for(int i = 0; i< cutLhs.size(); i++) {
				addLocal(ilcplex.ge(ilcplex.diff(cutLhs.get(i),cutRhs.get(i)), 0));
			}
		}
	}
	

	public static void main(String[] args) throws IOException, IloException {
		double U1=5;
		double U2=3;
		int num1=5;
		int num2=10;
		
		double[][] w= {{3,2,5,7,2},{6,2,4,9,3},{3,1,6,5,2},{4,1,3,2,2},{2,8,3,6,2},{7,5,4,9,3},{3,7,5,2,6},{6,4,2,4,6},{4,6,2,8,1},{3,4,8,5,3}};
		
		buildModel(w,num1,num2,U1,U2);
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161

下面为代码的求解结果：

参考文献

[1] Duran M A , Grossmann I E . An outer-approximation algorithm for a class of mixed-integer nonlinear programs[J]. Mathematical Programming, 1986, 36(3):307-339.