外点法惩罚函数（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();
}
}
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);
}*/
}
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;
}
}
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();
}*/
}
原文地址：http://blog.csdn.net/i_love_home/article/details/36242529
