i<sqrt(N)+1和i<=sqrt(N)不能相互替换?

博客通过实例展示了在C++中,使用for循环与sqrt函数时,`i<=sqrt(N)`与`i<sqrt(N)+1`的微妙差异。由于sqrt函数返回double类型,可能导致结果不一致。文章通过代码演示并分析了这种现象,强调了精度问题对编程结果的影响。

在我们平常的认知中,初学者一般都会认为在for循环中,i<N+1和i<=N是可以相互替换的,所以也会认为i<sqrt(N)+1和i<=sqrt(N)也是可以相互替换的,下面我来举几个例子来说明他们是否可以相互替换。

#include<iostream>
#include<cmath>
using namespace std;
int main()
{
    int N,count=0,num=0;
    cin>>N;
    for(int i=0;i<=sqrt(N);i++)
    {
    	count+=i;
	}
	for(int i=0;i<sqrt(N)+1;i++)
    {
    	num+=i;
	}
	cout<<count<<endl<<num;
    return 0;
}

我们先对N赋值一个10:

结果显示如此小的数它们就不相等了,我们继续分别把i值和开方值打印出来看一下

代码段:

#include<iostream>
#include<cmath>
using namespace std;
int main()
{
    int N,count=0,num=0;
    cin>>N;
    for(int i=0;i<=sqrt(N);i++)
    {
    	cout<<"i="<<i<<"|sqrt(N)="<<sqrt(N)<<" ";
    	count+=i;
	}
	cout<<endl;
	for(int i=0;i<sqrt(N)+1;i++)
    {
    	cout<<"i="<<i<<"|sqrt(N)="<<sqrt(N)<<" ";
    	num+=i;
	}
    return 0;
}

结果段:

 发现一个 一个执行到了i=3,一个执行到了i=4;上部分:3<=3.16228而下部分:4<3.16228+1;

说明sqrt函数返回的是double类型,就会导致题中不能替换的结果。

sqrt源码函数:

 

 

module constants implicit none integer, parameter :: dp = 8 ! 双精度 integer, parameter :: nion = 5 ! 3 Na? + 2 Cl? integer, parameter :: ndim = 3 * nion ! 总自由度 integer, parameter :: nfree = ndim - 6 ! 自由离子自由度 end module constants module utils use constants implicit none contains function vector_dot(a, b, n) result(dot) integer, intent(in) :: n real(dp), intent(in) :: a(n), b(n) real(dp) :: dot integer :: i dot = 0.0d0 do i = 1, n dot = dot + a(i)*b(i) end do end function vector_dot subroutine matrix_vector_product(A, x, y, n) integer, intent(in) :: n real(dp), intent(in) :: A(n,n), x(n) real(dp), intent(out) :: y(n) integer :: i, j y = 0.0d0 do i = 1, n do j = 1, n y(i) = y(i) + A(i,j)*x(j) end do end do end subroutine matrix_vector_product subroutine outer_product(a, b, C, n) integer, intent(in) :: n real(dp), intent(in) :: a(n), b(n) real(dp), intent(out) :: C(n,n) integer :: i, j do i = 1, n do j = 1, n C(i,j) = a(i) * b(j) end do end do end subroutine outer_product end module utils module chemistry use constants use utils implicit none contains subroutine unpack_coordinates(x, rfix, r) real(dp), intent(in) :: x(nfree), rfix(6) real(dp), intent(out) :: r(ndim) real(dp) :: d0 d0 = 2.5d0 ! 固定 Na1 Na2 r(1:3) = rfix(1:3) r(4:6) = rfix(4:6) ! Na3: 构成等边三角形 r(7) = d0 * 0.5d0 r(8) = d0 * sqrt(3.0d0) * 0.5d0 r(9) = 0.0d0 ! Cl1 Cl2 r(10:12) = x(1:3) r(13:15) = x(4:6) end subroutine unpack_coordinates subroutine total_energy(r, energy) real(dp), intent(in) :: r(ndim) real(dp), intent(out) :: energy real(dp) :: dr(3), rij, vij integer :: charge(nion) integer :: i, j data charge /1, 1, 1, -1, -1/ ! Na+, Na+, Na+, Cl-, Cl- energy = 0.0d0 do i = 1, nion - 1 do j = i + 1, nion dr(1) = r(3*i - 2) - r(3*j - 2) dr(2) = r(3*i - 1) - r(3*j - 1) dr(3) = r(3*i ) - r(3*j ) rij = sqrt(sum(dr**2)) if (rij < 1.0d-8) then energy = 1.0d10 return end if if (charge(i) * charge(j) < 0) then vij = -1.09d3 * exp(-rij / 0.321d0) / rij else vij = 1.09d3 * exp(-rij / 0.321d0) / rij endif vij = vij + 1.0d0 * (0.1d0 / rij)**12 energy = energy + vij end do end do end subroutine total_energy subroutine set_fixed_coords(rfix) real(dp), intent(out) :: rfix(6) real(dp) :: d0 d0 = 2.5d0 rfix = (/ 0.0d0, 0.0d0, 0.0d0, d0, 0.0d0, 0.0d0 /) end subroutine set_fixed_coords end module chemistry module optimization use constants use utils use chemistry implicit none contains subroutine numerical_gradient(x, rfix, grad) real(dp), intent(in) :: x(nfree), rfix(6) real(dp), intent(out) :: grad(nfree) real(dp), parameter :: h = 1.0d-8 real(dp) :: x_plus(nfree), x_minus(nfree), r(ndim), f_plus, f_minus integer :: i do i = 1, nfree x_plus = x; x_plus(i) = x(i) + h x_minus = x; x_minus(i) = x(i) - h call unpack_coordinates(x_plus, rfix, r) call total_energy(r, f_plus) call unpack_coordinates(x_minus, rfix, r) call total_energy(r, f_minus) grad(i) = (f_plus - f_minus) / (2.0d0 * h) end do end subroutine numerical_gradient subroutine update_hessian(H, s, y, gamma, n) integer, intent(in) :: n real(dp), intent(inout) :: H(n,n) real(dp), intent(in) :: s(n), y(n), gamma real(dp) :: Hs(n), temp1(n,n), temp2(n,n) call matrix_vector_product(H, s, Hs, n) call outer_product(s, s, temp1, n) temp1 = temp1 * (1.0d0 + vector_dot(y, Hs, n)/gamma) / gamma call outer_product(Hs, Hs, temp2, n) temp2 = temp2 / vector_dot(y, Hs, n) H = H + temp1 - temp2 end subroutine update_hessian subroutine bfgs_update(x, g, H, n) integer, intent(in) :: n real(dp), intent(inout) :: x(n), g(n), H(n,n) real(dp), parameter :: alpha_max = 1.0d0, c1 = 1.0d-6, rho = 0.5d0 real(dp) :: f0, fp, x_try(n), r(ndim), rfix(6) real(dp) :: alpha, s(n), y(n), gamma, g_new(n) integer :: iter call set_fixed_coords(rfix) call unpack_coordinates(x, rfix, r) call total_energy(r, f0) ! 回溯线搜索 alpha = alpha_max do while (alpha > 1.0d-10) x_try = x - alpha * g call unpack_coordinates(x_try, rfix, r) call total_energy(r, fp) if (fp <= f0 - c1 * alpha * vector_dot(g, g, n)) then exit else alpha = alpha * rho end if end do if (alpha <= 1.0d-10) then return end if s = -alpha * g x = x + s call numerical_gradient(x, rfix, g_new) y = g_new - g gamma = vector_dot(s, y, n) if (gamma > 1.0d-10) then call update_hessian(H, s, y, gamma, n) g = g_new end if end subroutine bfgs_update end module optimization program main use constants use chemistry use optimization implicit none integer, parameter :: maxiter = 200 real(dp), parameter :: eps = 1.0d-5 real(dp) :: x(nfree), g(nfree), H(nfree,nfree) real(dp) :: r(ndim), rfix(6), energy integer :: iter, i ! ?? 显式声明 i logical :: converged ! 初始化 call set_fixed_coords(rfix) call set_initial_guess(x) call numerical_gradient(x, rfix, g) ! 初始化 Hessian:单位阵 × 0.1 H = 0.0d0 do i = 1, nfree H(i,i) = 0.1d0 end do ! 主循环 converged = .false. do iter = 1, maxiter call unpack_coordinates(x, rfix, r) call total_energy(r, energy) write(*,'(A,I3,A,F15.8,A,E10.3)') & 'Iter ', iter, ': E=', energy, ' |grad|=', sqrt(vector_dot(g, g, nfree)) call bfgs_update(x, g, H, nfree) if (vector_dot(g, g, nfree) < eps**2) then converged = .true. exit end if end do ! 输出结果 if (converged) then write(*,'(/,A,I0,A/)') '? Converged in ', iter, ' iterations.' else write(*,'(/,A/)') '? Optimization did not converge!' end if call unpack_coordinates(x, rfix, r) write(*,'(A)') 'Final coordinates:' do i = 1, nion write(*,'(A,I1,2X,F9.5,2X,F9.5,2X,F9.5)') & 'Ion', i, r(3*i-2), r(3*i-1), r(3*i) end do ! 导出 XYZ 文件 open(unit=10, file='na3cl2.xyz', status='replace') write(10,'(I0)') nion write(10,'(A)') 'Na3Cl2 optimized structure' do i = 1, nion if (i <= 3) then write(10,'(A,3F10.5)') 'Na', r(3*i-2), r(3*i-1), r(3*i) else write(10,'(A,3F10.5)') 'Cl', r(3*i-2), r(3*i-1), r(3*i) end if end do close(10) write(*,'(/,A)') 'XYZ file saved: na3cl2.xyz' pause contains subroutine set_initial_guess(x) real(dp), intent(out) :: x(nfree) real(dp) :: d0, zcl d0 = 2.5d0; zcl = 1.0d0 x(1:3) = (/ d0/2, d0*sqrt(3.0d0)/6, zcl /) x(4:6) = (/ d0/2, d0*sqrt(3.0d0)/6, -zcl /) end subroutine set_initial_guess end program main,再代码里添加输出Na_Cl键长
11-02
tic; clc; clear; %% 海底管道多相流(油-水-气)仿真 - 多工况版本(含保温层) - 优化版 % 主要改进:引入漂移通量模型、压力校正法、强化稳定性控制 % ==================== 工程参数设置 ==================== % 仿真控制参数 t_end = 10*3600; % 仿真总时间 (s) dt = 0.1; % 初始时间步长 (s),后续动态调整 max_dt = 10; % 最大时间步长 (s) % 工况选择 (1-稳定生产, 2-停输, 3-再启动) condition = 1; % 工况编号 % 边界条件 T_in = 323; % 入口温度 (K) - 50°C Pin = 10e6; % 入口压力 (Pa) - 10 MPa Pout = 8e6; % 出口压力 (Pa) - 8 MPa [重新引入出口压力边界] u_in = 2; % 入口流速 (m/s) % 管道参数 L = 50e3; % 管道长度 (m) D = 0.3; % 管道内径 (m) epsilon = 0.0001; % 管道粗糙度 (m) N = 200; % 空间离散点数 CFL = 0.5; % CFL数 % 保温层参数 insulation_enabled = 0; % 是否启用保温层 (1-启用, 0-不启用) insulation_thickness = 0.05; % 保温层厚度 (m) k_insulation = 0.05; % 保温层导热系数 (W/m·K) % 各相物性参数 rho_o = 850; % 油相密度 (kg/m3) rho_w = 1025; % 水相密度 (kg/m3) rho_g = 200; % 气相密度 (kg/m3) mu_o = 0.05; % 油相粘度 (Pa·s) mu_w = 0.001; % 水相粘度 (Pa·s) mu_g = 1.2e-5; % 气相粘度 (Pa·s) alpha_o_in = 0.6; % 入口油相体积分数 alpha_w_in = 0.3; % 入口水相体积分数 alpha_g_in = 0.1; % 入口气相体积分数 % 热力学参数 T_sea = 277; % 海水温度 (K) - 4°C U_sea = 5; % 海水传热系数 (W/m2·K) c_o = 2000; % 油相比热容 (J/kg·K) c_w = 4180; % 水相比热容 (J/kg·K) c_g = 2200; % 气相比热容 (J/kg·K) k_pipe = 50; % 管道导热系数 (W/m·K) k_cond = 0.15; % 流体导热系数 (W/m·K) % 漂移通量模型参数 (关键改进:考虑气相滑移) C0 = 1.2; % 分布系数 Vgj = 0.3; % 漂移速度 (m/s) [针对气-液系统] % ==================== 参数设置结束 ==================== % 初始化空间网格 dx = L/(N-1); x = linspace(0, L, N); % 计算总传热系数 (修正热阻计算) if insulation_enabled == 1 D_insulation = D + 2*insulation_thickness; % 保温层外径 R_sea = 1 / (U_sea * pi * D_insulation); % 海水侧热阻 (K/W/m) R_insulation = log(D_insulation/D) / (2 * pi * k_insulation); R_total = R_sea + R_insulation; U_total = 1 / (R_total * pi * D); % 折算到管道内径的总传热系数 else U_total = U_sea; end % 初始化变量 u = u_in*ones(1, N); % 混合流速(m/s) alpha_o = alpha_o_in * ones(1, N); alpha_w = alpha_w_in * ones(1, N); alpha_g = alpha_g_in * ones(1, N); P = linspace(Pin, Pout, N); % 初始压力线性分布 T = T_in * ones(1, N); % 创建图形窗口 figure('Position', [100, 100, 1200, 800]); subplot(3, 2, 1); plot_o = plot(x, alpha_o, 'b', 'DisplayName', '油相'); hold on; plot_w = plot(x, alpha_w, 'g', 'DisplayName', '水相'); plot_g = plot(x, alpha_g, 'r', 'DisplayName', '气相'); xlabel('管道位置 (m)'); ylabel('体积分数'); title('各相体积分数随管道位置变化'); legend('Location', 'southoutside', 'Orientation', 'horizontal'); grid on; hold off; subplot(3, 2, 2); plot_P = plot(x, P/1e5, 'k'); xlabel('管道位置 (m)'); ylabel('压力 (bar)'); title('压力随管道位置变化'); grid on; subplot(3, 2, 3); plot_u = plot(x, u, 'y'); xlabel('管道位置 (m)'); ylabel('流速 (m/s)'); title('流速随管道位置变化'); grid on; subplot(3, 2, 4); plot_T = plot(x, T-273.15, 'r'); % 显示摄氏度 xlabel('管道位置 (m)'); ylabel('温度 (°C)'); title('温度分布'); grid on; subplot(3, 2, 6); time_text = text(0.1, 0.9, '时间: 0 s', 'Parent', gca, 'Units', 'normalized'); axis off; % 主循环 last_output_time = 0; output_interval = max(t_end/100, dt*10); % 输出间隔 % 用于压力校正的变量 P_corr = zeros(1, N); % 压力校正量 divergence = zeros(1, N); % 速度散度 for t = 0:dt:t_end % 动态调整时间步长 (强化CFL条件) u_max = max(abs(u)) + 1e-6; c_sound = sqrt(1.4 * max(P) / min([rho_o, rho_w, rho_g])); % 声速估计 dt_conv = CFL * dx / (u_max + c_sound + 1e-6); % 粘性扩散时间步长限制 mu_max = max([mu_o, mu_w, mu_g]); rho_min = min([rho_o, rho_w, rho_g]); dt_visc = 0.25 * dx^2 * rho_min / (mu_max + 1e-6); % 热扩散时间步长限制 cp_min = min([c_o, c_w, c_g]); k_eff = k_cond + 4*U_total*dx; % 有效导热系数 dt_therm = 0.25 * dx^2 * rho_min * cp_min / (k_eff + 1e-6); % 选择最严格的时间步长 dt = min([dt_conv, dt_visc, dt_therm, max_dt]); dt = max(dt, 1e-6); % 避免负时间步长 % 再启动工况下流速控制 if condition == 3 if t < t_end/5 u_in_current = min(u_in, u_in * t / (t_end/5)); else u_in_current = u_in; end else u_in_current = u_in; end % === 计算混合物性质 === rho_mix = alpha_o.*rho_o + alpha_w.*rho_w + alpha_g.*rho_g; rho_mix = max(rho_mix, 1e-6); % 计算混合物比热容 cp_mix = (alpha_o.*rho_o*c_o + alpha_w.*rho_w*c_w + alpha_g.*rho_g*c_g) ./ rho_mix; cp_mix = max(cp_mix, 1e-6); % 粘度模型 mu_mix = mu_o.*alpha_o + mu_w.*alpha_w + mu_g.*alpha_g; mu_mix = max(mu_mix, 1e-6); % 雷诺数计算 Re = rho_mix .* abs(u) * D ./ mu_mix; Re = max(Re, 1e-6); % 摩擦系数 (使用Haaland近似公式,避免迭代) f = zeros(size(Re)); for i = 1:length(Re) if Re(i) < 2300 f(i) = 64 / Re(i); % 层流 else % Haaland近似公式 f(i) = (-1.8 * log10( (epsilon/D/3.7)^1.11 + 6.9/Re(i) ) )^(-2); f(i) = max(f(i), 0.001); % 最小值限制 end end % === 相分数方程 (关键改进:引入漂移通量模型) === alpha_o_old = alpha_o; alpha_w_old = alpha_w; alpha_g_old = alpha_g; % 液相总体积分数 (油+水) alpha_l = 1 - alpha_g; alpha_l = max(alpha_l, 1e-6); % 计算气相漂移速度 v_g = C0 * u + Vgj; % 气相速度 for i = 2:N-1 % 油相(迎风格式) if u(i) >= 0 flux_o = u(i) * alpha_o_old(i-1); else flux_o = u(i) * alpha_o_old(i); end alpha_o(i) = alpha_o_old(i) - dt/dx * (flux_o - (u(i-1)>=0)*u(i-1)*alpha_o_old(i-1) + (u(i-1)<0)*u(i-1)*alpha_o_old(i)); % 水相(迎风格式) if u(i) >= 0 flux_w = u(i) * alpha_w_old(i-1); else flux_w = u(i) * alpha_w_old(i); end alpha_w(i) = alpha_w_old(i) - dt/dx * (flux_w - (u(i-1)>=0)*u(i-1)*alpha_w_old(i-1) + (u(i-1)<0)*u(i-1)*alpha_w_old(i)); % 气相(漂移通量模型) if v_g(i) >= 0 flux_g = v_g(i) * alpha_g_old(i-1); else flux_g = v_g(i) * alpha_g_old(i); end alpha_g(i) = alpha_g_old(i) - dt/dx * (flux_g - (v_g(i-1)>=0)*v_g(i-1)*alpha_g_old(i-1) + (v_g(i-1)<0)*v_g(i-1)*alpha_g_old(i)); end % 边界条件 alpha_o(1) = alpha_o_in; alpha_w(1) = alpha_w_in; alpha_g(1) = alpha_g_in; % 相分数限制归一化 alpha_o = max(min(alpha_o, 1.0), 0); alpha_w = max(min(alpha_w, 1.0), 0); alpha_g = max(min(alpha_g, 1.0), 0); alpha_sum = alpha_o + alpha_w + alpha_g; alpha_sum = max(alpha_sum, 1e-6); alpha_o = alpha_o ./ alpha_sum; alpha_w = alpha_w ./ alpha_sum; alpha_g = alpha_g ./ alpha_sum; % 更新混合物密度 rho_mix = alpha_o.*rho_o + alpha_w.*rho_w + alpha_g.*rho_g; rho_mix = max(rho_mix, 1e-6); % === 动量方程 (显式求解) === u_old = u; for i = 2:N-1 % 压力梯度项 (中心差分) dPdx = (P(i-1) - P(i+1))/(2*dx); % 摩擦项 friction = -0.5 * f(i) * rho_mix(i) * u_old(i)*abs(u_old(i)) / D; % 重力项 (假设管道水平) gravity = 0; % 对流项 (迎风) if u_old(i) >= 0 conv_term = -u_old(i) * (u_old(i) - u_old(i-1))/dx; else conv_term = -u_old(i) * (u_old(i+1) - u_old(i))/dx; end % 动量方程 dudt = conv_term + (dPdx + friction + gravity) / rho_mix(i); u(i) = u_old(i) + dt * dudt; end % 边界条件 u(1) = u_in_current; % 入口流速 u(N) = u(N-1); % 出口外推 % 限制流速范围 u = max(min(u, 100), -100); % === 压力校正 (SIMPLE算法简化版) === % 计算速度散度 for i = 2:N-1 divergence(i) = (u(i+1) - u(i-1)) / (2*dx); end % 压力泊松方程 (简化迭代) P_old = P; for iter = 1:3 % 简单迭代3次 for i = 2:N-1 P_corr(i) = 0.5 * (P_corr(i-1) + P_corr(i+1) - dx^2 * rho_mix(i) * divergence(i) / dt); end end % 更新压力速度 P(2:N-1) = P_old(2:N-1) + P_corr(2:N-1); % 速度校正 for i = 2:N-1 u(i) = u(i) - dt/rho_mix(i) * (P_corr(i+1) - P_corr(i-1))/(2*dx); end % 压力边界条件 P(1) = Pin; P(N) = Pout; % 压力限制 P = max(min(P, 100e6), 1e3); % === 能量方程 === T_old = T; for i = 2:N-1 % 热对流项 (迎风格式) if u(i) >= 0 dTdx = (T_old(i) - T_old(i-1)) / dx; else dTdx = (T_old(i+1) - T_old(i)) / dx; end % 热传导项 (中心差分) d2Tdx2 = (T_old(i-1) - 2*T_old(i) + T_old(i+1)) / dx^2; % 管壁热交换 q_wall = (4*U_total/D) * (T_sea - T_old(i)); % 能量方程 dTdt = -u(i)*dTdx + (k_cond/(rho_mix(i)*cp_mix(i)))*d2Tdx2 + q_wall/(rho_mix(i)*cp_mix(i)); T(i) = T_old(i) + dt * dTdt; end % 温度边界条件 if condition == 1 T(1) = T_in; % 稳定生产时入口温度保持恒定 else T(1) = T_sea; % 停输再启动时入口温度为环境温度 end T(N) = T(N-1); % 出口绝热边界 % 确保温度合理 T = max(min(T, 1000), 100); % === 更新图形输出 === if t - last_output_time >= output_interval || t == 0 || t >= t_end % 更新曲线 set(plot_o, 'YData', alpha_o); set(plot_w, 'YData', alpha_w); set(plot_g, 'YData', alpha_g); set(plot_P, 'YData', P/1e5); set(plot_u, 'YData', u); set(plot_T, 'YData', T-273.15); set(time_text, 'String', ['时间: ', num2str(t, '%.0f'), ' s']); % 根据工况输出 switch condition case 1 fprintf('稳定生产 - 时间: %.0f s, 入口压力: %.2f bar, 出口温度: %.2f °C\n', ... t, P(1)/1e5, T(end)-273.15); case 2 fprintf('停输 - 时间: %.0f s, 平均温度: %.2f °C\n', ... t, mean(T)-273.15); case 3 fprintf('再启动 - 时间: %.0f s, 入口流速: %.2f m/s, 出口压力: %.2f bar\n', ... t, u_in_current, P(end)/1e5); end drawnow; last_output_time = t; end end toc; %% 结果分析 fprintf('\n=== 仿真结果分析 ===\n'); fprintf('工况: '); switch condition case 1 fprintf('稳定生产\n'); case 2 fprintf('停输\n'); case 3 fprintf('再启动\n'); end if insulation_enabled == 1 fprintf('保温层状态: 启用\n'); else fprintf('保温层状态: 未启用\n'); end fprintf('总传热系数: %.3f W/m²·K\n', U_total); fprintf('入口压力: %.2f bar\n', P(1)/1e5); fprintf('出口压力: %.2f bar\n', P(end)/1e5); fprintf('压力降: %.2f bar\n', (P(1)-P(end))/1e5); fprintf('入口温度: %.2f °C\n', T(1)-273.15); fprintf('出口温度: %.2f °C\n', T(end)-273.15); fprintf('温度降: %.2f °C\n', T(1)-T(end)); fprintf('平均流速: %.2f m/s\n', mean(abs(u))); 该代码里流速压力不收敛,且超出正常工程范围
08-14
VectorVec4d SpatiotemporalQPOptimizer::optimize(const VectorVec4d& init_path) { size_t N = init_path.size(); if (N < 3) { std::cerr << "[QP WARNING] Path too short, skip optimization.\n"; return init_path; } const int vars_per_point = 4; const int total_vars = N * vars_per_point; // 1. 构造 Hessian Gradient Eigen::SparseMatrix<double> hessian(total_vars, total_vars); Eigen::VectorXd gradient = Eigen::VectorXd::Zero(total_vars); std::vector<Eigen::Triplet<double>> hessian_triplets; double w_pos = 10.0, w_theta = 1.0, w_time = 1.0; auto add_symmetric_triplet = [&](int row, int col, double val) { if (row >= 0 && row < total_vars && col >= 0 && col < total_vars) { hessian_triplets.emplace_back(row, col, val); if (row != col) hessian_triplets.emplace_back(col, row, val); // 对称项 } }; // Add small diagonal terms for ALL variables (regularization) for (size_t i = 0; i < N; ++i) { int idx = i * vars_per_point; add_symmetric_triplet(idx + 0, idx + 0, 1e-6); // x add_symmetric_triplet(idx + 1, idx + 1, 1e-6); // y add_symmetric_triplet(idx + 2, idx + 2, 1e-6); // theta add_symmetric_triplet(idx + 3, idx + 3, 1e-6); // t } // Add smoothing terms for intermediate points for (size_t i = 1; i < N - 1; ++i) { int idx = i * vars_per_point; // Position smoothing (x, y) for (int d = 0; d < 2; ++d) { int center = idx + d; int prev = center - vars_per_point; int next = center + vars_per_point; add_symmetric_triplet(center, center, 2 * w_pos); add_symmetric_triplet(prev, center, -w_pos); add_symmetric_triplet(next, center, -w_pos); } // Orientation smoothing (theta) int center_theta = idx + 2; int prev_theta = center_theta - vars_per_point; int next_theta = center_theta + vars_per_point; add_symmetric_triplet(center_theta, center_theta, 2 * w_theta); add_symmetric_triplet(prev_theta, center_theta, -w_theta); add_symmetric_triplet(next_theta, center_theta, -w_theta); // Time smoothing (t) int center_time = idx + 3; int prev_time = center_time - vars_per_point; int next_time = center_time + vars_per_point; add_symmetric_triplet(center_time, center_time, 2 * w_time); add_symmetric_triplet(prev_time, center_time, -w_time); add_symmetric_triplet(next_time, center_time, -w_time); } // 构造稀疏矩阵(合并重复 triplet) hessian.setFromTriplets(hessian_triplets.begin(), hessian_triplets.end(), [](double a, double b) { return a + b; }); hessian.makeCompressed(); // Debug: Print Hessian dimensions std::cout << "Hessian size: " << hessian.rows() << " x " << hessian.cols() << std::endl; // 校验 if (hessian.rows() != total_vars || hessian.cols() != total_vars) { std::cerr << "[QP ERROR] Hessian size mismatch.\n"; return init_path; } Eigen::MatrixXd H_dense(hessian); double asym_norm = (H_dense - H_dense.transpose()).norm(); if (asym_norm > 1e-6) { std::cerr << "[QP ERROR] Hessian is not symmetric! Norm: " << asym_norm << std::endl; return init_path; } if (hessian.nonZeros() == 0) { std::cerr << "[QP WARN] Hessian is zero! Applying small regularization.\n"; hessian.insert(0, 0) = 1e-6; hessian.makeCompressed(); } // 2. 构造约束 Eigen::SparseMatrix<double> constraint_mat; Eigen::VectorXd lower_bound, upper_bound; std::vector<Eigen::Triplet<double>> constraint_triplets; std::vector<double> lower, upper; for (size_t i = 0; i < N; ++i) { int base = i * vars_per_point; // 时间递增约束:t[i] > t[i-1] + 0.01 constraint_triplets.emplace_back(lower.size(), base + 3, 1.0); lower.push_back(i == 0 ? init_path[i][3] : init_path[i - 1][3] + 0.01); upper.push_back(init_path[i][3] + 1.0); // 可调节最大步长 // 方向变化限制:|theta[i] - theta[i-1]| <= max_theta_rate_ if (i > 0) { constraint_triplets.emplace_back(lower.size(), base + 2, 1.0); constraint_triplets.emplace_back(lower.size(), base + 2 - vars_per_point, -1.0); lower.push_back(-max_theta_rate_); upper.push_back(max_theta_rate_); } // 碰撞提示(但不约束) if (isInCollision(init_path[i][0], init_path[i][1])) { std::cerr << "[WARN] Init path point " << i << " in obstacle.\n"; } } constraint_mat.resize(lower.size(), total_vars); constraint_mat.setFromTriplets(constraint_triplets.begin(), constraint_triplets.end()); constraint_mat.makeCompressed(); lower_bound = Eigen::Map<Eigen::VectorXd>(lower.data(), lower.size()); upper_bound = Eigen::Map<Eigen::VectorXd>(upper.data(), upper.size()); // 3. 设置 OSQP OsqpEigen::Solver solver; solver.settings()->setVerbosity(false); solver.settings()->setWarmStart(true); solver.settings()->setPolish(true); solver.data()->setNumberOfVariables(total_vars); solver.data()->setNumberOfConstraints(lower_bound.size()); if (!solver.data()->setHessianMatrix(hessian)) { std::cerr << "[QP ERROR] Failed to set Hessian matrix.\n"; return init_path; } if (!solver.data()->setGradient(gradient)) { std::cerr << "[QP ERROR] Failed to set gradient.\n"; return init_path; } if (!solver.data()->setLinearConstraintsMatrix(constraint_mat)) { std::cerr << "[QP ERROR] Failed to set constraint matrix.\n"; return init_path; } if (!solver.data()->setLowerBound(lower_bound)) { std::cerr << "[QP ERROR] Failed to set lower bounds.\n"; return init_path; } if (!solver.data()->setUpperBound(upper_bound)) { std::cerr << "[QP ERROR] Failed to set upper bounds.\n"; return init_path; } if (!solver.initSolver()) { std::cerr << "[QP ERROR] Failed to initialize solver.\n"; return init_path; } auto result = solver.solveProblem(); if (result != OsqpEigen::ErrorExitFlag::NoError) { std::cerr << "[QP ERROR] Solver failed. Code: " << static_cast<int>(result) << std::endl; return init_path; } Eigen::VectorXd QPSolution = solver.getSolution(); if (QPSolution.size() != total_vars) { std::cerr << "[QP ERROR] Solution size mismatch.\n"; return init_path; } // 4. 构造输出路径 VectorVec4d optimized_path; for (size_t i = 0; i < N; ++i) { int idx = i * vars_per_point; Eigen::Vector4d pt; pt << QPSolution[idx], QPSolution[idx + 1], QPSolution[idx + 2], QPSolution[idx + 3]; optimized_path.push_back(pt); } std::cerr << "[INFO] QP optimization complete.\n"; return optimized_path; } 修改这个函数里面的代码,解决如下问题Hessian size: 1544 x 1544 ERROR in LDL_factor: Error in KKT matrix LDL factorization when computing the nonzero elements. The problem seems to be non-convex ERROR in osqp_setup: KKT matrix factorization. The problem seems to be non-convex. [OsqpEigen::Solver::initSolver] Unable to setup the workspace. [QP ERROR] Failed to initialize solver.
08-05
%% Sobol敏感性分析实现 - 修正版本 % 处理58行10列数据集,前7列为输入变量,后3列为输出变量 % 使用正确的Sobol敏感性分析方法 clear; clc; close all; %% 主函数 fprintf('Sobol敏感性分析程序启动...\n'); fprintf('==================================================\n'); try % 1. 加载数据 data = load_user_data(); fprintf('\n数据加载完成: %d行 × %d列\n', size(data,1), size(data,2)); fprintf('数据预览:\n'); % 创建变量名 varNames = [strseq('X',1:7), strseq('Y',1:3)]; dataTable = array2table(data, 'VariableNames', varNames); disp(head(dataTable)); % 2. 准备分析数据 [X, Y, problem] = prepare_sobol_analysis(data); % 3. 执行真正的Sobol分析 results = perform_true_sobol_analysis(X, Y, problem); % 4. 识别高敏感性变量 high_sensitivity_vars = identify_high_sensitivity_variables(results, 0.1); % 5. 可视化结果 visualize_sobol_results(results); fprintf('\n==================================================\n'); fprintf('分析完成!\n'); fprintf('==================================================\n'); catch ME fprintf('分析过程中出现错误: %s\n', ME.message); rethrow(ME); end %% 辅助函数:生成序列名称 function names = strseq(prefix, indices) names = arrayfun(@(i) sprintf('%s%d', prefix, i), indices, 'UniformOutput', false); end %% 加载用户数据 function data = load_user_data() fprintf('正在加载数据集...\n'); try data = readmatrix('分析数据.csv'); fprintf('成功从 筛选后的数据.xlsx 加载数据\n'); catch fprintf('无法加载数据文件,请检查文件是否存在或格式是否正确\n'); fprintf('将使用示例数据进行演示\n'); data = generate_sample_data(); end % 验证数据维度 if size(data, 2) ~= 10 fprintf('警告: 数据列数不是10,当前为%d列\n', size(data, 2)); fprintf('请确保数据格式为:前7列输入,后3列输出\n'); end % 检查数据中是否有NaN或Inf值 if any(isnan(data(:)) | isinf(data(:))) fprintf('警告: 数据中包含NaN或Inf值,将尝试处理...\n'); % 用列均值替换NaN for i = 1:size(data, 2) col_data = data(:, i); nan_indices = isnan(col_data); if any(nan_indices) col_mean = mean(col_data(~nan_indices & ~isinf(col_data))); data(nan_indices, i) = col_mean; fprintf(' 列 %d: 替换了 %d 个NaN值\n', i, sum(nan_indices)); end % 处理Inf值 inf_indices = isinf(col_data); if any(inf_indices) col_median = median(col_data(~inf_indices & ~isnan(col_data))); data(inf_indices, i) = col_median; fprintf(' 列 %d: 替换了 %d 个Inf值\n', i, sum(inf_indices)); end end end end %% 生成示例数据(仅在无法加载用户数据时使用) function data = generate_sample_data() fprintf('使用示例数据进行演示...\n'); % 生成示例数据集:58行10列,前7列输入,后3列输出 rng(42); % 设置随机种子 % 生成输入变量 (X1-X7) n_samples = 58; n_inputs = 7; % 创建输入变量 - 使用均匀分布 X = rand(n_samples, n_inputs); % 生成输出变量 (Y1-Y3),基于输入变量的非线性组合,包含交互效应 Y1 = 2.0 * X(:, 1) + 1.5 * X(:, 2).^2 + 0.5 * X(:, 3) .* X(:, 4) + 0.1 * randn(n_samples, 1); Y2 = 1.0 * X(:, 2) + 3.0 * X(:, 5).^3 + 0.8 * X(:, 6) + 0.2 * X(:, 1) .* X(:, 7) + 0.1 * randn(n_samples, 1); Y3 = 0.5 * X(:, 1) .* X(:, 7) + 2.0 * X(:, 2) + 1.0 * X(:, 6).^2 + 0.3 * X(:, 3) .* X(:, 5) + 0.1 * randn(n_samples, 1); % 组合数据 data = [X, Y1, Y2, Y3]; end %% 准备Sobol分析 function [X, Y, problem] = prepare_sobol_analysis(data) % 提取输入输出 X = data(:, 1:7); % 前7列是输入 Y = data(:, 8:10); % 后3列是输出 % 检查输出变量是否有零方差 for i = 1:size(Y, 2) if var(Y(:, i)) == 0 fprintf('警告: 输出变量 Y%d 方差为零,添加微小噪声...\n', i); Y(:, i) = Y(:, i) + 1e-10 * randn(size(Y, 1), 1); end end % 计算每个输入变量的实际范围 bounds = [min(X); max(X)]'; % 定义问题 problem = struct(); problem.num_vars = 7; problem.names = strseq('X', 1:7); problem.bounds = bounds; % 使用实际数据的范围 fprintf('输入变量范围:\n'); for i = 1:7 fprintf(' %s: [%.4f, %.4f]\n', problem.names{i}, bounds(i,1), bounds(i,2)); end end %% 执行真正的Sobol分析 function results = perform_true_sobol_analysis(X, Y, problem) fprintf('\n正在执行真正的Sobol敏感性分析...\n'); results = struct(); % 对每个输出变量进行分析 for i = 1:size(Y, 2) output_name = sprintf('Y%d', i); fprintf('\n分析输出变量: %s\n', output_name); y = Y(:, i); % 使用基于方差的Sobol方法 try % 方法1: 使用Saltelli方法近似计算Sobol指数 [S1, ST] = calculate_sobol_indices(X, y); catch ME fprintf('Sobol分析失败: %s\n', ME.message); fprintf('使用替代方法...\n'); % 方法2: 使用基于回归的方法 [S1, ST] = calculate_regression_based_indices(X, y); end % 创建结果结构 results.(output_name) = struct(); results.(output_name).S1 = S1; % 一阶敏感性指数 results.(output_name).ST = ST; % 总效应指数 results.(output_name).names = problem.names; % 打印结果 fprintf('\n%s的Sobol敏感性分析结果:\n', output_name); fprintf('变量\t\t一阶指数\t总效应指数\t差值(交互效应)\n'); for j = 1:length(problem.names) interaction_effect = ST(j) - S1(j); fprintf('%s\t\t%.4f\t\t%.4f\t\t%.4f\n', problem.names{j}, S1(j), ST(j), interaction_effect); end % 计算并显示交互效应占比 total_interaction = sum(ST - S1); fprintf('总交互效应: %.4f\n', total_interaction); end end %% 计算Sobol指数 - 基于Saltelli方法 function [S1, ST] = calculate_sobol_indices(X, y) n_vars = size(X, 2); n_samples = size(X, 1); % 初始化结果 S1 = zeros(n_vars, 1); ST = zeros(n_vars, 1); % 计算总方差 total_variance = var(y); if total_variance == 0 error('输出变量方差为零,无法计算Sobol指数'); end % 对每个变量计算一阶效应总效应 for i = 1:n_vars % 一阶效应:固定其他变量,只变化第i个变量 % 使用条件方差估计 % 分组数据:按第i个变量的值分组 [sorted_X, sort_idx] = sort(X(:, i)); sorted_y = y(sort_idx); % 使用移动窗口计算条件期望 window_size = max(3, floor(n_samples / 10)); % 窗口大小 conditional_means = zeros(n_samples, 1); for j = 1:n_samples start_idx = max(1, j - floor(window_size/2)); end_idx = min(n_samples, j + floor(window_size/2)); conditional_means(j) = mean(sorted_y(start_idx:end_idx)); end % 计算一阶效应方差 variance_S1 = var(conditional_means); S1(i) = variance_S1 / total_variance; % 总效应:使用基于抽样的方法 % 固定第i个变量,随机重排其他变量 n_permutations = min(100, n_samples); total_effect_variances = zeros(n_permutations, 1); for p = 1:n_permutations % 创建新样本:固定第i个变量,重排其他变量 X_permuted = X; other_vars = setdiff(1:n_vars, i); % 重排其他变量 perm_idx = randperm(n_samples); X_permuted(:, other_vars) = X(perm_idx, other_vars); % 使用KNN预测输出 if n_samples > 10 k = min(5, floor(n_samples/10)); y_pred = knn_predict(X, y, X_permuted, k); total_effect_variances(p) = var(y_pred); else % 样本太少,使用简单方法 total_effect_variances(p) = total_variance; end end % 计算总效应 ST(i) = mean(total_effect_variances) / total_variance; end % 确保指数在合理范围内 S1 = max(0, min(1, S1)); ST = max(S1, min(1, ST)); % 总效应应该不小于一阶效应 end %% 使用KNN进行预测 function y_pred = knn_predict(X_train, y_train, X_test, k) n_test = size(X_test, 1); y_pred = zeros(n_test, 1); for i = 1:n_test % 计算距离 distances = sqrt(sum((X_train - X_test(i, :)).^2, 2)); % 找到最近的k个邻居 [~, idx] = sort(distances); nearest_indices = idx(1:min(k, length(idx))); % 使用邻居的平均值作为预测 y_pred(i) = mean(y_train(nearest_indices)); end end %% 基于回归的方法计算敏感性指数 function [S1, ST] = calculate_regression_based_indices(X, y) n_vars = size(X, 2); % 标准化数据 X_scaled = zscore(X); y_scaled = zscore(y); % 拟合线性模型获取一阶效应 mdl = fitlm(X_scaled, y_scaled); coefficients = mdl.Coefficients.Estimate(2:end); % 排除截距 % 一阶效应:标准化回归系数的平方 S1 = (coefficients .^ 2) / sum(coefficients .^ 2); % 总效应:使用随机森林重要性作为近似 try % 使用随机森林估计总效应 tree = fitrensemble(X_scaled, y_scaled, 'Method', 'Bag', 'NumLearningCycles', 50); imp = oobPermutedPredictorImportance(tree); ST = imp / sum(imp); catch % 如果随机森林失败,使用线性模型加上交互项 ST = calculate_total_effect_with_interactions(X_scaled, y_scaled); end % 确保总效应不小于一阶效应 ST = max(S1, ST); end %% 计算包含交互效应的总效应 function ST = calculate_total_effect_with_interactions(X, y) n_vars = size(X, 2); ST = zeros(n_vars, 1); for i = 1:n_vars % 创建包含所有变量与第i个变量交互项的模型 interaction_terms = []; for j = 1:n_vars if j ~= i interaction_terms = [interaction_terms, X(:, i) .* X(:, j)]; end end % 组合特征 features = [X, interaction_terms]; % 拟合模型 mdl = fitlm(features, y); % 计算第i个变量及其交互项的贡献 var_contrib = 0; for j = 1:length(mdl.CoefficientNames) coef_name = mdl.CoefficientNames{j}; if contains(coef_name, sprintf('x%d', i)) || ... (contains(coef_name, 'x') && contains(coef_name, sprintf('x%d_', i))) var_contrib = var_contrib + mdl.Coefficients.Estimate(j)^2; end end ST(i) = var_contrib; end % 归一化 ST = ST / sum(ST); end %% 识别高敏感性变量 function high_sensitivity_vars = identify_high_sensitivity_variables(results, threshold) if nargin < 2 threshold = 0.1; end fprintf('\n%s\n', repmat('=', 1, 60)); fprintf('高敏感性变量识别结果:\n'); fprintf('%s\n', repmat('=', 1, 60)); high_sensitivity_vars = struct(); output_names = fieldnames(results); for i = 1:length(output_names) output = output_names{i}; result = results.(output); fprintf('\n%s的高敏感性变量 (阈值: %.1f):\n', output, threshold); % 基于总效应指数识别高敏感性变量 high_indices = find(result.ST > threshold & ~isnan(result.ST)); if ~isempty(high_indices) high_sensitivity_vars.(output) = {}; for j = 1:length(high_indices) idx = high_indices(j); var_name = result.names{idx}; S1 = result.S1(idx); ST = result.ST(idx); interaction = ST - S1; fprintf(' %s: 一阶=%.4f, 总效应=%.4f, 交互效应=%.4f\n', ... var_name, S1, ST, interaction); high_sensitivity_vars.(output){end+1} = {var_name, S1, ST, interaction}; end else fprintf(' 无变量超过阈值\n'); end end end %% 可视化Sobol结果 function visualize_sobol_results(results) fprintf('\n正在生成Sobol敏感性分析图表...\n'); output_names = fieldnames(results); n_outputs = length(output_names); % 创建主图 fig = figure('Position', [100, 100, 1400, 1000]); for i = 1:n_outputs output = output_names{i}; result = results.(output); subplot(2, 2, i); % 准备数据 variables = result.names; s1_values = result.S1; st_values = result.ST; interaction_values = st_values - s1_values; % 检查是否有NaN值 valid_indices = ~isnan(s1_values) & ~isnan(st_values); if ~all(valid_indices) fprintf('警告: %s 中有NaN值,将跳过这些变量\n', output); variables = variables(valid_indices); s1_values = s1_values(valid_indices); st_values = st_values(valid_indices); interaction_values = interaction_values(valid_indices); end if isempty(variables) text(0.5, 0.5, '无有效数据', 'HorizontalAlignment', 'center', 'FontSize', 14); title(sprintf('%s - 无有效数据', output), 'FontSize', 12); continue; end % 创建堆叠柱状图显示一阶效应交互效应 x = 1:length(variables); % 绘制一阶效应 bars1 = bar(x, s1_values, 0.8, 'FaceColor', [0.2, 0.6, 0.8]); hold on; % 绘制交互效应(堆叠在一阶效应之上) bars2 = bar(x, interaction_values, 0.8, 'FaceColor', [0.8, 0.4, 0.2]); % 设置堆叠 bars2.BaseValue = 0; % 添加总效应线 plot(x, st_values, 'ko-', 'LineWidth', 2, 'MarkerSize', 6, 'MarkerFaceColor', 'k'); xlabel('输入变量'); ylabel('敏感性指数'); title(sprintf('%s的Sobol敏感性分析', output), 'FontSize', 12); set(gca, 'XTick', x, 'XTickLabel', variables); xtickangle(45); % legend({'一阶效应 (S1)', '交互效应', '总效应 (ST)'}, 'Location', 'best'); grid on; % 添加阈值线 yline(0.1, '--r', '阈值=0.1', 'LabelHorizontalAlignment', 'left', 'FontSize', 8); % 添加数值标签 for j = 1:length(x) text(x(j), st_values(j) + 0.02, sprintf('ST=%.3f', st_values(j)), ... 'HorizontalAlignment', 'center', 'VerticalAlignment', 'bottom', 'FontSize', 8); end end % 添加综合比较图在第四个位置 if n_outputs >= 3 subplot(2, 2, 4); % 准备比较数据 s1_matrix = zeros(n_outputs, length(results.(output_names{1}).names)); st_matrix = zeros(n_outputs, length(results.(output_names{1}).names)); for i = 1:n_outputs s1_matrix(i, :) = results.(output_names{i}).S1; st_matrix(i, :) = results.(output_names{i}).ST; end % 检查是否有NaN值 if any(isnan(s1_matrix(:))) || any(isnan(st_matrix(:))) fprintf('警告: 比较图数据中包含NaN值,将用0替换\n'); s1_matrix(isnan(s1_matrix)) = 0; st_matrix(isnan(st_matrix)) = 0; end % 创建分组柱状图 x = 1:size(s1_matrix, 2); width = 0.35; for i = 1:n_outputs bar_positions = x + (i-1) * width - (n_outputs-1)*width/2; bar(bar_positions, st_matrix(i, :), width, 'FaceAlpha', 0.7); hold on; end xlabel('输入变量'); ylabel('总敏感度指数 (ST)'); title('各输出变量的输入变量总敏感性比较', 'FontSize', 12); set(gca, 'XTick', x, 'XTickLabel', results.(output_names{1}).names); % legend(output_names, 'Location', 'best'); grid on; end % 调整布局 sgtitle('Sobol敏感性分析结果 (S1: 一阶效应, ST: 总效应)', 'FontSize', 16, 'FontWeight', 'bold'); % 保存图表 saveas(fig, 'sobol_sensitivity_analysis_corrected.png'); fprintf('图表已保存为: sobol_sensitivity_analysis_corrected.png\n'); end上面的代码运行之后得到的结果我认为存在问题,得到的一阶灵敏度指标与全局灵敏度指标相差值过大,且得到的结果准确性我认为不强,请帮我对代码进行修改完善,保证运行之后得到的结果更准确
最新发布
11-22
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

小浪浪、

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值