C++实现坐标轮换法；黄金分割一维搜索；外推内插初始区间(最优化计算)

编程实现坐标轮换法，子问题求解采用外推内插+抛物线逼近 或 外推内插+黄金分割。

外推内插法用来确定初始搜索区间。

黄金分割用来缩小搜索区间，最终取区间中心作为一维搜索结果。

main.cpp

~main.cpp
void run();
void waitAndExit();int main() {run();waitAndExit();
}

imp.cpp

#include<iostream>
using std::cout;
using std::endl;#include<vector>
using std::vector;#include<cassert>#include<string>
using std::string;#include<iomanip>const float cc = 0.01f; //cc for convergence condition
const float cc_gold = 0.01f;
const float golde = 0.618034f;
const int dim = 3;float target(float x1, float x2, float x3) {return 3 * x1 * x1 + 2 * x2 * x2 + 1 * x3 * x3;
}float target(vector<float>x) {return target(x[0], x[1], x[2]);
}float target(vector<float>x, float dimdata, int dim) {x[dim] = dimdata;return target(x);
}float norm2(vector<float>x) {float result = 0.0f;for (int i = 0; i < x.size(); i++) {result += x[i] * x[i];}return sqrtf(result);
}vector<float> substract(vector<float> x, vector<float> y) {//perform r = x -y;vector<float> r;size_t length = x.size();r.resize(length);for (int i = 0; i < length; i++) {r[i] = x[i] - y[i];}return r;
}void printV(vector<float>x, string prefix = "",int indent = 0) {while (indent-->0){cout << "\t";}if (prefix.size() > 0)cout << prefix << " ";for (int i = 0; i < x.size(); i++) {cout << x[i] << " ";}cout << endl;
}void swap(float&x, float&y) {float tmp = x;x = y;y = tmp;
}void run() {cout.setf(std::ios::fixed);cout << std::setprecision(6); //set output formatsvector<float> X0(dim, 0);vector<float> Xn(dim, 0);vector<float> Dx(dim, 1e9f); // to hold Xn+1 -Xnfor (int i = 0; i < dim; i++) X0[i] = i + 1; //given X0, See 8-1, Page 119Xn = X0;// to hold the result from Extrapolation and interpolation method, "left" and "right" to initialize golden split intervalfloat left;	 float right;	float center; int current_dim = 0; //for cyclic coordinate method, update by dimx = (dimx +1 )% totaldimdo {cout << "DIM[" << current_dim+1 << "]" << endl;float step = 1;float factor = 1;float x1 = Xn[current_dim];float x2 = x1 + step * factor++;float x3 = x2 + step * factor++;float xkprevprev = x1;float xkprev = x2;float xk = x3;//extrapolation------------------------------------------------------------------cout << "\textra~ and inter~ method:" << endl;cout << "\tinit " << xkprevprev << " " << xkprev << " " << xk << endl;if (target(Xn, xkprevprev, current_dim) > target(Xn, xkprev, current_dim)) {//do nothing}else if (target(Xn, xkprevprev, current_dim) < target(Xn, xkprev, current_dim)) {step = -1; factor = 1; //reset factorxkprev = xkprevprev + step * factor++;xk = xkprev + step * factor++;}else {//unkonwn error}while (target(Xn, xk, current_dim) < target(Xn, xkprev, current_dim)) {xkprevprev = xkprev;xkprev = xk;xk = xk + step * factor++;cout << "\tproc " << xkprevprev << "\t" << xkprev << "\t" << xk << endl;}//extrapolation------------------------------------------------------------------// interpolation------------------------------------------------------------------float xknext = (xk + xkprev) / 2.f;if (target(Xn, xknext, current_dim) < target(Xn, xkprev, current_dim)) {left = xkprev;center = xknext;right = xk;}else if (target(Xn, xknext, current_dim) > target(Xn, xkprev, current_dim)){left = xkprevprev;center = xkprev;right = xknext;}else {//unkonwn error}if (left > right) swap(left, right);cout << "\tinte " << left << "\t" << center<< "\t" << right << endl;// interpolation------------------------------------------------------------------//golden split------------------------------------------------------------------cout << "\n\t" << "golden split" << endl;while (right - left > cc_gold) {float r_split = left + golde * (right - left);float l_split = left + (1 - golde) * (right - left); //need to modify herecout << "\t" << left << "\t" << l_split << "\t" << r_split << "\t" << right << endl;if (target(Xn, l_split, current_dim) < target(Xn, r_split, current_dim)) {right = r_split;}else{left = l_split;}}float xi = (left + right) / 2.f;  //golden split result//golden split------------------------------------------------------------------vector<float>Xnnext = Xn;Xnnext[current_dim] = xi; //update Xn+1Dx = substract(Xnnext, Xn); // compute DxprintV(Xn, "Xn");printV(Xnnext, "Xn+1");cout << "F(X) = " << target(Xnnext) << "        <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<" <<  endl;printV(Dx, "DX");cout << "Norm2DX " << norm2(Dx) << endl << endl << endl;Xn = Xnnext;current_dim = (current_dim + 1) % dim; //switch to next dim } while (norm2(Dx) > cc);printV(Xn, "Final");cout << "F(X) = " << target(Xn) << endl;return;
}void waitAndExit() {cout << "Press Enter to exit." << endl;getchar();return;
}

cslayee#163.com