fatique_curve.cpp - main - starlv/Fatigue_plan_optimizer_cpp

fatique_curve.cpp
198 строк6 KB
Agamirov L.V.
Add files via upload
22 май 2026, 18:40
22 май 2026, 18:401f16d79
#include <cmath>
#include <vector>
#include <iostream>
#include <fstream>
#include <iomanip>
#include <map> 
#include <numeric>
#include <algorithm>
#include <functional>

using namespace std;

int neldermead(vector<double>& x0, double eps, double(*func)(vector<double>));
int fatique_curve_fit(int ku,vector<double> cx,vector<int> ni,vector<double> lgN,vector<double> slgN,
                       double &sigma_inf,double &C,double &m,double &Q);


struct fatique{
 vector<double>w;
 vector<double>cx;
 vector<double> lgN;
 int ku;
};
fatique fat;


//###################################################################

int neldermead(vector<double>&x0, double eps,double(*func)(vector<double>)) {

    double rho,chi, psi, sigma, nonzdelt, zdelt;
    double maxfun,fval,fxr,fxe,fxc,fxcc;
    int i,j,k,N,maxiter,iterations,doshrink;

    rho = 1;chi = 2;psi = 0.5;sigma = 0.5;nonzdelt = 0.05;zdelt = 0.00025;
    N = x0.size();maxiter = N * 300; maxfun = 1.e20;

    vector<vector<double>> sim(N + 1,vector<double>(N,0.0));
    vector<double> fsim(N + 1);
    sim[0]= x0;
    vector<double>y = x0;

    for (k = 0; k < N; k++) {
        if (y[k] != 0) {
            y[k] = (1 + nonzdelt) * y[k];
        } else {
            y[k] = zdelt;
        }
        sim[k + 1] = y;
    }
  
    for (i = 0; i < N + 1; i++)   fsim[i] = func(sim[i]);

    iterations = 0;
    fval = maxfun;

//###########################################################

    while (true) {
        if (fval <= eps || iterations >= maxiter) break;
        vector<double> xbar(N,0.0);
        for (i = 0; i < N; i++) {
            for ( j = 0; j < N; j++) xbar[j] += sim[i][j];
        }
        for (i = 0; i < N; i++) xbar[i] /= N;
        
        vector<double> xr(N);
        for (i = 0; i < N; i++) xr[i] = (1 + rho) * xbar[i] - rho * sim[N][i];
        fxr = func(xr);

        doshrink = 0;
        if (fxr < fsim[0]) {
            vector<double> xe(N);
            for (i = 0; i < N; i++)  xe[i] = (1 + rho * chi) * xbar[i] - rho * chi * sim[N][i];
            fxe = func(xe);
            if (fxe < fxr) {
                sim[N] = xe;fsim[N] = fxe;
            } else {
                sim[N] = xr;fsim[N] = fxr;
            }
        } else {
            if (fxr < fsim[N - 1]) {
                sim[N] = xr;fsim[N] = fxr;
            } else {
                if (fxr < fsim[N]) {
                    vector<double> xc(N);
                    for (i = 0; i < N; i++) xc[i] = (1 + psi * rho) * xbar[i] - psi * rho * sim[N][i];
                    fxc = func(xc);
                    if (fxc <= fxr) {
                        sim[N] = xc; fsim[N] = fxc;
                    } else {
                        doshrink = 1;
                    }
                } else {
                    vector<double> xcc(N);
                    for (i = 0; i < N; i++)  xcc[i] = (1 - psi) * xbar[i] + psi * sim[N][i];
                    fxcc = func(xcc);
                    if (fxcc < fsim[N]) {
                        sim[N] = xcc;fsim[N] = fxcc;
                    } else {
                        doshrink = 1;
                    }
                }
                if (doshrink) {
                    for (j = 1; j < N + 1; j++) {
                        for (i = 0; i < N; i++)     sim[j][i] = sim[0][i] + sigma * (sim[j][i] - sim[0][i]);
                        fsim[j] = func(sim[j]);
                    }
                }
            }
        }

        iterations++;
        fval = func(sim[0]);
        vector<int> ind(N + 1);
        for (i = 0; i < N + 1; i++) ind[i] = i;
        sort(ind.begin(), ind.end(), [&](int i, int j) { return fsim[i] < fsim[j]; });
        vector<vector<double>> sim_sorted(N + 1, vector<double>(N));
        vector<double> fsim_sorted(N + 1);
        for (i = 0; i < N + 1; i++) {
            sim_sorted[i] = sim[ind[i]];
            fsim_sorted[i] = fsim[ind[i]];
        }
        sim = sim_sorted; fsim = fsim_sorted;

    }  //end while

    x0=sim[0];
    return iterations;

}

//===========Целевая функция==============================================


double target_fun(vector<double> params_norm) {

   int i,k;
   double sigma_inf,C,m,sum_sq,z,sigma_adj,pred,residual;

    k=fat.cx.size();
    sigma_inf=params_norm[0];
    C=params_norm[1];
     m=params_norm[2];
    
    // Границы (в нормализованных координатах)
    if (sigma_inf<0) return 1e20;
    if (sigma_inf>=fat.cx[k-1]) return 1e20;
    if (m<=0.0) return 1e20;
    if (C<=0.0) return 1e20;

    
    sum_sq=0.0;z=0;
    for (i=0;i<k; ++i) {
        pred=(log10(C)-log10(fat.cx[i]-sigma_inf))/m;
        if(fat.ku==0) residual=pred-log10(fat.lgN[i]);
        if(fat.ku==1) residual=(pred-fat.lgN[i]);
        sum_sq+=fat.w[i]*residual*residual;
        z+=fat.w[i];
    }
    return sum_sq/z;
}


//==============================================================

int fatique_curve_fit(int ku,vector<double> cx,vector<int> ni,vector<double> lgN,vector<double> slgN,
                       double &sigma_inf,double &C,double &m,double &Q) {

    int i,k,iterations;
    double z,eps; 

    //ku=0; Cx=sigma_inf+C*(lgN)^(-m) 
    //ku=1; Cx=sigma_inf+C*(N)^(-m)


    k=cx.size();
    vector<double> w;
    fat.cx=cx;fat.lgN=lgN;fat.ku=ku;
    vector<double> params_norm = {sigma_inf,C,m};
    eps=1e-4;
    for (i=0;i<k;++i) {
          if(ku==0) z=ni[i]*lgN[i]*lgN[i]/(slgN[i]*slgN[i]);
          if(ku==1) z=ni[i];///(slgN[i]*slgN[i]);
          w.push_back(z);
    }
    fat.w=w;

    iterations=neldermead(params_norm,eps,target_fun);
    
    Q=target_fun(params_norm);
    sigma_inf=params_norm[0];
    C=params_norm[1];
    m=params_norm[2];

    return iterations;
  
}