外点惩罚函数法·约束优化问题

外点法惩罚函数（r增加，SUMT.java）用于求解约束优化问题，解题步骤如下：

Step1 输入目标函数与约束方程，构建外点惩罚函数法求解方程，求解初始化。

Step2 对求解方程进行一次无约束优化方法求解（鲍威尔BWE），得到新解。

Step3 新解与原解求误差，如误差满足精度要求，则输出解，否则增加因子r，执行Step 2。

鲍威尔法（BWE.java）是N维无约束求解方法，需要调用一维求解方法，一维求解方法采用黄金分割法（GSM.java）。

在实现算法的代码中，我去掉了输入处理，人为地将输入确定下来，可减少代码篇幅。

我会将文件打包放入我的下载，欢迎大家一起交流。

（1）外点法惩罚函数 SUMT.java：

package ODM.Method;

import java.util.Arrays;
/*
 * 无约束优化方法：惩罚函数法·外点法
 */
public class SUMT {

	private int n = 6;							// 维数，变量个数
	private final double eps = 1e-5;			// 精度
	private final double c = 5;				// 递增系数
	private double r = 0.1;					// 惩罚因子，趋向无穷
	public SUMT(){
		Finit();
		AlgorithmProcess();
		AnswerOutput();
	}
	// 结果
	private double[] xs;
	private double ans;
	private void Finit(){
		xs = new double[n];
		Arrays.fill(xs, 0);
		ans = -1;
		//xs[0] = xs[1] = xs[2] = xs[4] = 1;	xs[3] = 3;	xs[5] = 5;
	}
	// 算法主要流程
	private void AlgorithmProcess(){
		int icnt = 0;						// 迭代次数
		double[] x = new double[n];		// 转化为无约束优化问题的解
		while(true){
			icnt++;
			BWE temp = new BWE(n, r, xs);	// 采用鲍威尔方法求函数最优解
			x = temp.retAns();
			if(retOK(x) <= eps){			// 满足精度要求
				for(int i = 0; i < n; i++)
					xs[i] = x[i];
				ans = temp.mAns();
				break;
			}
			r = c * r;
			for(int i = 0; i < n; i++)
				xs[i] = x[i];
		}
		System.out.println("迭代次数：" + icnt);
	}
	// 收敛条件（只有一个，不完善）
	private double retOK(double[] x){
		double sum = 0;
		for(int i = 0; i < n; i++){
			sum += Math.pow(x[i] - xs[i], 2);
		}
		return Math.sqrt(sum);
	}
	// 结果输出
	private void AnswerOutput(){
		for(int i = 0; i < n; i++)
			System.out.printf("%.6f\t", xs[i]);
		System.out.printf("%.6f\n", ans);
	}

	public static void main(String[] args) {
		// TODO Auto-generated method stub
		new SUMT();
	}

}

（2）鲍威尔法 BWE.java：

package ODM.Method;

import java.util.Arrays;

public class BWE {

	private double r;
	// 初始化变量
	private double[] x0;						// 初始解集
	private double[][] e;						// 初始方向
	private int N;
	final private double eps = 1e-5;
	private Func F;

	// 初始化：初始点, 初始矢量(n 个，n*n 矩阵), 维数
	private void Init(int n){
		this.x0 = new double[n];
		if(r == -1)
			Arrays.fill(this.x0, 0);
		else{

		}

		this.e = new double[n][n];
		for(int i = 0; i < n; i++){
			for(int j = 0; j < n; j++){
				if(i != j)e[i][j] = 0;
				else e[i][j] = 1;
			}
		}

		this.N = n;
		if(r != -1)
			F = new Func(r);
		else
			F = new Func();
	}

	// 搜索点, 方向矢量
	private double[][] x;
	private double[][] d;

	// 方向重排, 队列操作
	private void queueDir(double[] X){
		// 删去首方向
		for(int i = 0; i < N-1; i++){
			for(int j = 0; j < N; j++){
				d[i][j] = d[i+1][j];
			}
		}
		// 新方向插入队尾
		for(int i = 0; i < N; i++)
			d[N-1][i] = X[i];
	}

	private void Process(){

		x = new double[N+1][N];
		d = new double[N][N];

		for(int j = 0; j < N; j++)
			x[0][j] = x0[j];
		for(int i = 0; i < N; i++){
			for(int j = 0; j < N; j++){
				d[i][j] = e[i][j];
			}
		}

		int k = 0;							// 迭代次数
		while(k < N){
			for(int i = 1; i <= N; i++){
				GSM t = new GSM(F, x[i-1], d[i-1]);
				x[i] = t.getOs();
			}
			double[] X = new double[N];
			for(int i = 0; i < N; i++)
				X[i] = x[N][i] - x[0][i];
			queueDir(X);
			GSM t = new GSM(F, x[N], X);
			x[0] = t.getOs();
			k++;
		}
	}
	// 答案打印
	private void AnswerOutput(){
		for(int i = 0; i < N; i++){
			System.out.printf("x[%d] = %.6f\n", i+1, x[0][i]);
//			System.out.print(x[0][i] + " ");
		}
		System.out.printf("最小值：%.6f\n", F.fGetVal(x[0]));
//		System.out.println(": " + F.fGetVal(x[0]));
	}

	public BWE(int n){
		this.r = -1;
		Init(n);
		Process();
		AnswerOutput();
	}

	public BWE(int n, double r, double[] x){
		this.r = r;
		Init(n);
		for(int i = 0; i < n; i++)
			x0[i] = x[i];
		Process();
	}
	// 返回结果，解向量和最优值
	public double[] retAns(){
		return x[0];
	}
	public double mAns(){
		return F.fGetVal(x[0], 0);
	}
	/*
	public static void main(String[] args) {
		// TODO Auto-generated method stub
		new BWE(2);
	}*/

}

（3）黄金分割 GSM.java：

package ODM.Method;
/*
 * 黄金分割法
 */
public class GSM {

	private int N;																// 维度
	private final double landa = (Math.sqrt(5)-1)/2;							// 0.618
	private double[] x1;
	private double[] x2;
	private double[] os;
	private final double eps = 1e-5;											// 解精度
	private ExtM EM;															// 用于获取外推法结果

	// 最优值输出
	public double[] getOs() {
		return os;
	}
	// 函数, 初始点, 方向矢量
	public GSM(Func Sample, double[] x, double[] e) {
		//for(int i = 0; i < e.length; i++)System.out.print(e[i] + " ");System.out.println();
		initial(Sample, x, e);
		process(Sample);
		AnswerPrint(Sample);
	}
	// 结果打印
	private void AnswerPrint(Func Sample) {
		os = new double[N];
		for(int i = 0; i < N; i++)
			os[i] = 0.5*(x1[i] + x2[i]);

//		System.out.println("os = " + os[0] + " " + os[1]);
//		System.out.println("ans = " + Sample.fGetVal(os));

	}
	// 向量范值
	private double FanZhi(double[] b, double[] a){
		double sum = 0;
		for(int i = 0; i < N; i++){
			if(b[i] - a[i] != 0 && b[i] == 0)return eps*(1e10);
			if(b[i] == 0)continue;
			sum += Math.pow((b[i] - a[i]) / b[i], 2);
		}
		return Math.pow(sum, 0.5);
	}
	// 算法主流程
	private void process(Func Sample) {
		double[] xx1 = new double[N];
		SubArraysCopy(xx1);
		double yy1 = Sample.fGetVal(xx1);

		double[] xx2 = new double[N];
		AddArraysCopy(xx2);
		double yy2 = Sample.fGetVal(xx2);
		// 迭代过程
		while(true){
			if(yy1 >= yy2){
				ArraysCopy(xx1, x1);
				ArraysCopy(xx2, xx1);	yy1 = yy2;
				AddArraysCopy(xx2);
				yy2 = Sample.fGetVal(xx2);
			}else{
				ArraysCopy(xx2, x2);
				ArraysCopy(xx1, xx2);	yy2 = yy1;
				SubArraysCopy(xx1);
				yy1 = Sample.fGetVal(xx1);
			}
			//System.out.println(FanZhi(x2, x1) + " / " + Math.abs((yy2 - yy1)/yy2));
			if(FanZhi(x2, x1) < eps && Math.abs(yy2 - yy1) < eps)break;
		}
	}
	// 获得外推法结果：左右边界
	private void initial(Func Sample, double[] x, double[] e) {
		N = x.length;
		EM = new ExtM(Sample, x, e);
		x1 = EM.getX1();
		x2 = EM.getX3();
	}
	// 向量赋值
	private void ArraysCopy(double[] s, double[] e){
		for(int i = 0; i < N; i++)
			e[i] = s[i];
	}
		// + landa
	private void AddArraysCopy(double[] arr){
		for(int i = 0; i < N; i++)
			arr[i] = x1[i] + landa*(x2[i] - x1[i]);
	}
		// - landa
	private void SubArraysCopy(double[] arr){
		for(int i = 0; i < N; i++)
			arr[i] = x2[i] - landa*(x2[i] - x1[i]);
	}
	/*
	public static void main(String[] args) {
		// TODO Auto-generated method stub
		double[] C = {0, 0};
		double[] d = {1, 0};
		new GSM(new Func(), C, d);
	}
	*/
}

以上算法文件包含函数方程，黄金分割时有一维搜索的外推法确定“高低高”区间。

函数方程 Func.java，外推法 ExtM.java。

Func.java：

package ODM.Method;

public class Func {

	private int N;								// N 维
	private double[] left;						// 函数左边界
	private double[] right;					// 函数右边界
	private double r;

	public Func() {
		r = -1;
	}
	public Func(double r) {
		this.r = r;
	}

	// 定义函数与函数值返回
	public double fGetVal(double[] x){
		if(r != -1)return fGetVal(x, r);
		// 10*(x1+x2-5)^2 + (x1-x2)^2
		return 10*Math.pow(x[0]+x[1]-5, 2) + Math.pow(x[0]-x[1], 2);
		//
	}
	private double max(double a, double b){return a > b ? a : b;}
	public double fGetVal(double[] x, double r){
		double ret = 0;
//		function f1
//		ret =  Math.pow(x[0]-5, 2) + 4*Math.pow(x[1]-6, 2)
//				+ r * (
//				+ Math.pow(max(64-x[0]*x[0]-x[1]*x[1], 0), 2)
//				+ Math.pow(max(x[1]-x[0]-10, 0),  0)
//				+ Math.pow(max(x[0]-10, 0), 2)
//				);

//		function f2
//		ret = x[0]*x[0] + x[1]*x[1] + r*(1-x[0]>0 ? 1-x[0] : 0)*(1-x[0]>0 ? 1-x[0] : 0);

//		function f3
		ret = Math.pow(x[0]-x[3], 2) + Math.pow(x[1]-x[4], 2) + Math.pow(x[2]-x[5], 2) +
				r * (
				+ Math.pow(max(x[0]*x[0]+x[1]*x[1]+x[2]*x[2]-5, 0), 2)
				+ Math.pow(max(Math.pow(x[3]-3, 2)+x[4]*x[4]-1, 0), 2)
				+ Math.pow(max(x[5]-8, 0), 2)
				+ Math.pow(max(4-x[5], 0), 2)
				);
		return ret;
	}
}

ExtM.java：

package ODM.Method;

/*
 * 外推法确定“高-低-高”区间
 */
public class ExtM {

	private int N;							// 函数维数
	private double[] x1;
	private double[] x2;
	private double[] x3;
	private double y1;
	private double y2;
	private double y3;
	private double h;						// 步长
	private double[] d;					// 方向矢量
	public double[] getX1() {
		return x1;
	}

	public double[] getX2() {
		return x2;
	}

	public double[] getX3() {
		return x3;
	}

	public double getH() {
		return h;
	}
	// 函数, 初始点，方向
	public ExtM(Func Sample, double[] x, double[] e) {
		initial(Sample, x, e);
		process(Sample);
		AnswerPrint();
	}
	// 初始化阶段
	private void initial(Func Sample, double[] x, double[] e)
	{
		N = x.length;
		x1 = new double[N];
		x2 = new double[N];
		x3 = new double[N];
		h = 0.01;
		d = new double[N];

		ArraysCopy(e, 0, d);
		//for(int i = 0; i < d.length; i++)System.out.print(d[i]);System.out.println();
		ArraysCopy(x, 0, x1);
		y1 = Sample.fGetVal(x1);

		ArraysCopy(x, h, x2);
		y2 = Sample.fGetVal(x2);

	}
	// 循环部分
	private void process(Func Sample)
	{
		if(y2 > y1){
			h = -h;
			ArraysCopy(x1, 0, x3);
			y3 = y1;
		}else{
			ArraysCopy(x2, h, x3);	y3 = Sample.fGetVal(x3);
		}
		while(y3 < y2){
			h = 2*h;
//			System.out.println("h = " + h);
			ArraysCopy(x2, 0, x1);		y1 = y2;
			ArraysCopy(x3, 0, x2);		y2 = y3;
			ArraysCopy(x2, h, x3);		y3 = Sample.fGetVal(x3);
//			System.out.println("x1 = " + x1[0] + " " + x1[1] + " y1 = " + y1);
//			System.out.println("x2 = " + x2[0] + " " + x2[1] + " y2 = " + y2);
//			System.out.println("x3 = " + x3[0] + " " + x3[1] + " y3 = " + y3);
		}

	}
	// 打印算法结果
	private void AnswerPrint()
	{
//		System.out.println("x1 = " + x1[0] + " " + x1[1] + " y1 = " + y1);
//		System.out.println("x2 = " + x2[0] + " " + x2[1] + " y2 = " + y2);
//		System.out.println("x3 = " + x3[0] + " " + x3[1] + " y3 = " + y3);
	}
	// 向量转移
	private void ArraysCopy(double[] s, double c, double[] e){
		for(int i = 0; i < s.length; i++)
			e[i] = d[i]*c + s[i];
	}

	/*
	public static void main(String[] args) {
		// TODO Auto-generated method stub
		// new ExtM();
	}*/

}

外点惩罚函数法·约束优化问题