工程优化学习(进退法、黄金分割法、二次插值法、三次插值法、最速下降法)

工程优化课程学了进退法、黄金分割法、二次插值法、三次插值法、最速下降法。我自己只测试用最速下降法求(恩，水货我只调试了这一个，调了好多次才调通，然后其他的也没心情调试了)：minf(x) = x1^2 + 25*x2^2，得到结果如下图所示，与该函数的极小值点（0，0）在精度误差之内。

#ifndef ENGINEERINGOPTIMIZATION_H_INCLUDED
#define ENGINEERINGOPTIMIZATION_H_INCLUDED

#include <vector>

using std::vector;

/* description:区间类
*/
class Region
{
public:
    double frontierLow; //一维定义域下界
    double frontierHigh; //一维定义域上界
    vector<double> lowerBoundDomain; //n维矢量下界
    vector<double> upperBoundDomain; //n维矢量上界

    Region() {}
    Region(double l, double h)
    {
        frontierLow = l;
        frontierHigh = h;
    }
    Region(double l, double h, vector<double> lowerBound,
            vector<double> upperBound)
    {
        frontierLow = l;
        frontierHigh = h;
        lowerBoundDomain = lowerBound;
        upperBoundDomain = upperBound;
    }
    Region(vector<double> lowerBound, vector<double> upperBound)
    {
        lowerBoundDomain = lowerBound;
        upperBound = upperBoundDomain;
    }
    Region(const Region &region)
    {
        frontierLow = region.frontierLow;
        frontierHigh = region.frontierHigh;
        lowerBoundDomain = region.lowerBoundDomain;
        upperBoundDomain = region.upperBoundDomain;
    }
};


/* @description 计算n维矢量的模值
**
*/
double vectorMag(const vector<double> &v);

/*
** @description 进退法求搜索区间。
** @param initialPoint 搜索初始点
** @param initialStep 搜索初始步长
** @param targetFunc 目标函数f(x)的回调函数，调用者需要提供此函数。
**  它接受一个参数，变量x，返回目标函数的函数值
** @return 返回搜索到的区间
*/
//变量x是一维空间中的点
Region computeSearchRegion(double initialPoint, double initialStep, double (*targetFuc)(double x, int derivationOrder));

/*
** @description 黄金分割法求最小值，目标函数需是单峰函数。缩短搜素区间的三个原则(a<=x<=b)：
**          1、“去坏留好”原则，若f(x1)<f(x2),则a<x_min<x2
**          2、对称原则，x1-a = b-x2
**          3、等比收缩原则，每次留下区间的长度都是上次留下区间长度的w倍(0<w<1)
**           调用者需提供目标函数
** @param region 目标函数的定义域
** @param precision 精度
** @param target 指向目标函数的函数指针
** @return 返回目标函数的极小值点
*/
double goldenSectionMethod(const Region &region, double precision, double (*targetFuc)(double x, int derivationOrder));


/* @description 二次插值法
** @return 返回目标函数的极小值点
*/
double quadraticInterpolation(double x1, double x2, double x3, double (*targetFuc)(double x, int derivationOrder));


/* @description 三次插值法（Davidon搜索法）
**              目标函数的第二参数为求导阶数
** @return 返回目标函数的极小值点
*/
double cubicInterpolation(const Region &region, double step, double precision,
                          double (*targetFuc)(double x, int derivationOrder));

/* @description 最速下降法求解极小值点 n维函数
** @return 返回极小值点n维矢量
** 回调函数 derivationOrder为负数时求最优步长，非负数求导
*/

vector<double> steepestDescent(vector<double> &initialPoint, double precision,
                                    vector<double> (*targetFuc)(vector<double> &x, int derivationOrder));


#endif // ENGINEERINGOPTIMIZATION_H_INCLUDED

#include <iostream>
#include "EngineeringOptimization.h"
#include <cassert>
#include <cmath>


/* @description 计算n维矢量的模值
**
*/
double vectorMag(const vector<double> &v)
{
    double mag = 0;
    for (int i = 0; i < v.size(); i++)
    {
        mag += v[i] * v[i];
    }

    return std::sqrt(mag);
}

/*
** @description 进退法求搜索区间。
** @param initialPoint 搜索初始点
** @param initialStep 搜索初始步长
** @param double (*targetFunc)(double, double) 目标函数f(x)的回调函数，调用者需要提供此函数。
**  它接受一个参数，变量x，返回目标函数的函数值
** @return 返回搜索到的区间
*/
Region computeSearchRegion(double initialPoint, double initialStep,
                           double (*targetFuc)(double x, int derivationOrder))
{
    double frontierLow = initialPoint;
    double frontierHigh = initialPoint + initialStep;
    double fucValueLow = targetFuc(frontierLow, 0);
    double fucValueHigh = targetFuc(frontierHigh, 0);
    double step = initialStep;
    bool FLAG = false;

    if (fucValueHigh < fucValueLow) //前进运算
    {
        do
        {
            if (FLAG)
                frontierLow -= step;

            step *= 2;
            frontierHigh += step;
            fucValueLow = fucValueHigh;
            fucValueHigh = targetFuc(frontierHigh, 0);
            FLAG = true;
        }
        while (fucValueHigh < fucValueLow);

        return Region(frontierLow, frontierHigh);
    }
    else //后退运算
    {
        step = -(step / 4);

        do
        {
            if (FLAG)
            {
                frontierHigh = frontierLow - step;
                step += step;
            }

            frontierLow += step;
            fucValueHigh = fucValueLow;
            fucValueLow = targetFuc(frontierLow, 0);
            FLAG = true;
        }
        while (fucValueHigh >= fucValueLow);

        return Region(frontierLow, frontierHigh);
    }

}


/*
** @description 黄金分割法求最小值，目标函数需是单峰函数。缩短搜素区间的三个原则(a<=x<=b)：
**          1、“去坏留好”原则，若f(x1)<f(x2),则a<x_min<x2
**          2、对称原则，x1-a = b-x2
**          3、等比收缩原则，每次留下区间的长度都是上次留下区间长度的w倍(0<w<1)
**           调用者需提供目标函数
** @param region 目标函数的定义域
** @param precision 精度
** @param target 指向目标函数的函数指针
** @return 返回目标函数的最小值
*/
double goldenSectionMethod(const Region &region, double precision,
                           double (*targetFuc)(double x, int derivationOrder))
{
    double frontierLow = region.frontierLow;
    double frontierHigh = region.frontierHigh;
    double x1 = frontierLow + 0.382 * (frontierHigh - frontierLow);
    double x2 = frontierLow + frontierHigh - x1;
    double valueLow = targetFuc(x1, 0);
    double valueHigh = targetFuc(x2, 0);

    while (valueLow > valueHigh)
    {
        frontierLow = x1;
        if (frontierHigh - frontierLow > -precision && frontierHigh - frontierLow < precision)
            return (frontierHigh + frontierLow) / 2.0;
        x1 = x2;
        x2 = frontierLow + 0.618 * (frontierHigh - frontierLow);
        valueLow = valueHigh;
        valueHigh = targetFuc(x2, 0);
    }

    while (valueLow <= valueHigh)
    {
        frontierHigh = x2;
        if (frontierHigh - frontierLow > -precision && frontierHigh - frontierLow < precision)
            return (frontierHigh + frontierLow) / 2.0;
        x2 = x1;
        x1 = frontierLow + 0.382 * (frontierHigh - frontierLow);
        valueHigh = valueLow;
        valueLow = targetFuc(x1, 0);
    }
}


/* @description 二次插值法
** @return 返回目标函数的最小值
*/
double quadraticInterpolation(double x1, double x2, double x3,
                              double (*targetFuc)(double x, int derivationOrder))
{
    double tmp;
    if (x1 < x2)
    {
        tmp = x1;
        x1 = x2;
        x2 = tmp;
    }
    if (x3 < x2)
    {
        tmp = x2;
        x2 = x3;
        x3 = tmp;
    }

    double f1 = targetFuc(x1, 0);
    double f2 = targetFuc(x2, 0);
    double f3 = targetFuc(x3, 0);

    double c1 = (f3 - f1) / (x3 - x1);
    double c2 = ((f2 - f1) / (x2 - x1) - c1) / (x2 - x3);
    double x_min = 0.5 * (x1 + x3 - c1 / c2);
    return x_min;
}


/* @description 三次插值法（Davidon搜索法）
**              目标函数的第二参数为求导阶数
** @return 返回目标函数最小值
*/
double cubicInterpolation(const Region &region, double step, double precision,
                          double (*targetFuc)(double x, int derivationOrder))
{
    double a = region.frontierLow;
    double b = region.frontierHigh;
    double h = step;
    bool FLAG = false;

    assert(precision > 0);

    double v = targetFuc(a, 1);
    if (v > -precision && v < precision)
        return a;

    do
    {
        if (FLAG)
            h /= 10;

        if (v < 0)
            h = std::abs(h);
        else
            h = -std::abs(h);

        double u;
        do
        {
            b = a + h;
            u = targetFuc(b, 1);
            if (u > - precision && u < precision)
                return b;
            if (u * v < 0)
                break;

            h *= 2;
            v =u;
            a = b;
        }
        while(true);

        double s = 3 * (targetFuc(a, 0) - targetFuc(b, 0)) / (b - a);
        double z = s - u -v;
        double w = std::sqrt(z * z - u * v);

        if (v > 0)
        {
            a = a - (b - a) * v / (z - w -v);
        }
        else
        {
            a = a - (b - a) * v / (z + w -v);
        }

        v = targetFuc(a, 1);

        FLAG = true;
    }
    while (v > precision || v < -precision);

    return a;
}


/* @description 最速下降法求解极小值点 n维函数
** @return 返回极小值点n维矢量
** 回调函数 derivationOrder为负数时求最优步长，非负数求导
*/
vector<double> steepestDescent(vector<double> &initialPoint, double precision,
                               vector<double> (*targetFuc)(vector<double> &x, int derivationOrder))
{
    vector<double> x_k;
    x_k = initialPoint;
    vector<double> x_k1;
    vector<double> s_k;
    double e = precision;
    do
    {
        s_k = targetFuc(x_k, 1);
        double mag = vectorMag(s_k);
        if (mag <= e)
            return x_k;

        vector<double> combine;
        for (int i = 0; i < x_k.size(); i++)
            combine.push_back(x_k.at(i));
        for (int i = 0; i < s_k.size(); i++)
            combine.push_back(s_k.at(i));

        vector<double> tmp = targetFuc(combine, -1);
        double step = vectorMag(tmp);

            std::cout << "step =" << step <<std::endl;

        x_k1.clear();
        for (int i = 0; i < x_k.size(); i++)
        {
            x_k1.push_back(x_k.at(i) - s_k.at(i) * step);
        }
        x_k = x_k1;

    } while(true);
}

#include <iostream>
#include <vector>
#include "EngineeringOptimization.h"

using namespace std;

vector<double> targetFuc(vector<double> &x, int derivationOrder);

int main()
{
    cout << "Hello world!" << endl;
    vector<double> x0;
    x0.push_back(2);
    x0.push_back(2);
    cout << "Intial Point: " << x0[0] << x0[1] <<endl;
    vector<double> min_value = steepestDescent(x0, 0.01, targetFuc);
    vector<double> result = targetFuc(min_value, 0);
    cout<< "x*= [" << min_value[0] <<" "<< min_value[1] <<"]"<<endl;;
    cout<<"The min value point is "<< result.at(0) <<endl;;
    return 0;
}

//min f(x)=x1^2 + 25x2^2
vector<double> targetFuc(vector<double> &x, int derivationOrder)
{
    vector<double> f;
    double value;
    f.push_back(0);
    switch (derivationOrder)
    {
    case 0:
        f.clear();
        value = x[0] * x[0] + 25 * x[1] * x[1];
        f.push_back(value);
        break;
    case 1:
        f.clear();
        value = 2 * x[0];
        f.push_back(value);
        value = 50 * x[1];
        f.push_back(value);
        break;
    default:
        break;
    }

    if (derivationOrder < 0)
    {
        f.clear();
        value = (2.0*x[0]*x[2]+50*x[1]*x[3])/(2*x[2]*x[2] + 50*x[3]*x[3]);
        f.push_back(value);
    }

    return f;

}