fatigue_curve.cpp - main - starlv/Fatigue_plan_optimizer_cpp

fatigue_curve.cpp
234 строки7 KB
Agamirov L.V.
Add files via upload
22 май 2026, 18:40
22 май 2026, 18:401f16d79
#include <unsupported/Eigen/NonLinearOptimization>
#include <Eigen/Core>
#include <Eigen/LU>
#include <vector>
#include <cmath>

using namespace std;

struct FatigueResult {
    double sigma_inf;
    double C;
    double m;
    double Q;
    Eigen::MatrixXd covariance;
    int optimization_status;
    int iterations;
};

struct FatigueFunctor {
    int m_inputs;
    int m_values;
    
    vector<double> cx;
    vector<double> lgN;
    vector<double> w;
    int ku;
    
    FatigueFunctor(int ku, const vector<double>& cx, 
                   const vector<double>& lgN, const vector<double>& w)
        : ku(ku), cx(cx), lgN(lgN), w(w) {
        m_inputs = 3;
        m_values = cx.size();
    }
    
    int operator()(const Eigen::VectorXd &x, Eigen::VectorXd &fvec) const {
        double sigma_inf = x(0);
        double C = pow(10.0, x(1));
        double m = x(2);
        
        fvec.resize(m_values);
        double pred, residual;
        for (int i = 0; i < m_values; ++i) {
            if (sigma_inf < 0 || sigma_inf >= cx[i] - 1.0 || m < 0.0 || m > 10.0) {
                fvec(i) = 1e10;
                continue;
            }
            
            pred = (log10(C) - log10(cx[i] - sigma_inf)) / m;
            if (ku == 0) residual = pred - log10(lgN[i]);
            if (ku == 1) residual = pred - lgN[i];
            fvec(i) = sqrt(w[i]) * residual;
        }
        return 0;
    }
    
    int df(const Eigen::VectorXd &x, Eigen::MatrixXd &fjac) const {
        const double eps = 1e-6;
        int n = m_values;
        int m = m_inputs;
        
        fjac.resize(n, m);
        
        Eigen::VectorXd f0(n);
        (*this)(x, f0);
        
        for (int j = 0; j < m; ++j) {
            Eigen::VectorXd x_plus = x;
            x_plus(j) += eps;
            
            Eigen::VectorXd f_plus(n);
            (*this)(x_plus, f_plus);
            
            for (int i = 0; i < n; ++i) {
                fjac(i, j) = (f_plus(i) - f0(i)) / eps;
            }
        }
        return 0;
    }
    
    void computeModelDerivatives(const Eigen::VectorXd &x, 
                                  Eigen::MatrixXd &J_model) const {
        double sigma_inf = x(0);
        double log10C = x(1);
        double m = x(2);
        
        J_model.resize(m_values, 3);
        
        for (int i = 0; i < m_values; ++i) {
            double sigma = cx[i];
            double delta = sigma - sigma_inf;
            
            if (delta <= 1e-10) {
                J_model.row(i).setZero();
                continue;
            }
            
            double pred = (log10C - log10(delta)) / m;
            
            // d(model)/d(sigma_inf) = 1/(m * ln(10) * delta)
            J_model(i, 0) = 1.0 / (m * log(10.0) * delta);
            
            // d(model)/d(log10C) = 1/m
            J_model(i, 1) = 1.0 / m;
            
            // d(model)/d(m) = -pred / m
            J_model(i, 2) = -pred / m;
        }
    }
    
    int inputs() const { return m_inputs; }
    int values() const { return m_values; }
};

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

vector<double> computeWeights(int ku,const vector<int>& ni,const vector<double>& lgN,const vector<double>& slgN) {

    size_t k = ni.size();
    vector<double> w(k);
    for (size_t i = 0; i < k; ++i) {
        if (ku == 0)
            w[i] = ni[i] * lgN[i] * lgN[i] / (slgN[i] * slgN[i]);
        else
            w[i] = ni[i] / (slgN[i] * slgN[i]);
    }
    return w;
}

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

void computeInitialParameters(int ku, const vector<double>& cx,const vector<double>& lgN,
                              double sigma_inf_init,double& logC_init, double& m_init) {

    double s1 = log10(cx[0] - sigma_inf_init);
    double s2 = log10(cx[1] - sigma_inf_init);
    
    if (ku == 0) {
        m_init = (s1 - s2) / (log10(lgN[1]) - log10(lgN[0]));
        logC_init = s1 + m_init * log10(lgN[0]);
    } else {
        m_init = (s1 - s2) / (lgN[1] - lgN[0]);
        logC_init = s1 + m_init * lgN[0];
    }
}

Eigen::MatrixXd computeCovarianceMatrix(FatigueFunctor& functor,
                                         const Eigen::VectorXd& x_final,
                                         double Q, int n_params) {
    // Вычисляем матрицу производных МОДЕЛИ
    Eigen::MatrixXd J_model;
    functor.computeModelDerivatives(x_final, J_model);
    
    // Строим диагональную матрицу весов
    int n_data = functor.m_values;
    Eigen::MatrixXd W = Eigen::MatrixXd::Zero(n_data, n_data);
    for (int i = 0; i < n_data; ++i) {
        W(i, i) = functor.w[i];
    }
    
    // Матрица нормальных уравнений: Jᵀ W J
    Eigen::MatrixXd H = J_model.transpose() * W * J_model;
    
    // Диагональная поправка для устойчивости
    double lambda = 1e-8;
    for (int i = 0; i < n_params; ++i)
        H(i, i) += lambda;
    
    // Обращение матрицы
    Eigen::MatrixXd covariance;
    Eigen::FullPivLU<Eigen::MatrixXd> lu(H);
    
    if (lu.isInvertible()) {
        covariance = lu.inverse();
        covariance *= Q;  // Масштабируем на дисперсию
    } else {
        covariance = Eigen::MatrixXd::Identity(n_params, n_params) * 1e10;
    }
    
    return covariance;
}

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

FatigueResult estimateFatigueCurve(int ku, const vector<double>& cx,
                                   const vector<int>& ni,
                                   const vector<double>& lgN,
                                   const vector<double>& slgN,
                                   double sigma_inf_guess) {
    FatigueResult result;
    result.optimization_status = -1;
    
    vector<double> w = computeWeights(ku, ni, lgN, slgN);
    
    double logC_init, m_init;
    computeInitialParameters(ku, cx, lgN, sigma_inf_guess, logC_init, m_init);
    
    FatigueFunctor functor(ku, cx, lgN, w);
    
    Eigen::VectorXd x(3);
    x << sigma_inf_guess, logC_init, m_init;
    
    Eigen::LevenbergMarquardt<FatigueFunctor, double> lm(functor);
    lm.parameters.maxfev = 2000;
    lm.parameters.ftol = 1e-12;
    lm.parameters.xtol = 1e-12;
    lm.parameters.gtol = 1e-12;
    
    int ret = lm.minimize(x);
    result.optimization_status = ret;
    result.iterations = lm.iter;
    
    result.sigma_inf = x(0);
    result.C = pow(10.0, x(1));
    result.m = x(2);
    
    // Вычисление Q
    double Q = 0.0;
    double z_total = 0.0;
    double pred, residual;
    
    for (size_t i = 0; i < cx.size(); ++i) {
        pred = (log10(result.C) - log10(cx[i] - result.sigma_inf)) / result.m;
        if (ku == 0) residual = pred - log10(lgN[i]);
        else residual = pred - lgN[i];
        Q += w[i] * residual * residual;
        z_total += w[i];
    }
    result.Q = Q / z_total;
    
    // Вычисляем ковариационную матрицу с использованием производных
    result.covariance = computeCovarianceMatrix(functor, x, result.Q, 3);
    
    return result;
}