%%% EKF UKF PF三种算法对比
clc
close all
clear;
% tic;
x = 0.1; % 初始状态
x_estimate = 1;%状态的估计
e_x_estimate = x_estimate; %EKF的初始估计
u_x_estimate = x_estimate; %UKF的初始估计
p_x_estimate = x_estimate; %PF的初始估计
Q = 10;%input('请输入过程噪声方差Q的值: '); % 过程状态协方差
R = 1;%input('请输入测量噪声方差R的值: '); % 测量噪声协方差
P =5;%初始估计方差
e_P = P; %UKF方差
u_P = P;%UKF方差
pf_P = P;%PF方差
tf = 50; % 模拟长度
x_array = [x];%真实值数组
e_x_estimate_array = [e_x_estimate];%EKF最优估计值数组
u_x_estimate_array = [u_x_estimate];%UKF最优估计值数组
p_x_estimate_array = [p_x_estimate];%PF最优估计值数组
u_k = 1; %微调参数
u_symmetry_number = 4; % 对称的点的个数
u_total_number = 2 * u_symmetry_number + 1; %总的采样点的个数
linear = 0.5;
N = 500; %粒子滤波的粒子数
close all;
%粒子滤波初始N 个粒子
for i = 1 : N
p_xpart(i) = p_x_estimate + sqrt(pf_P) * randn;
end
for k = 1 : tf
% 模拟系统
x = linear * x + (25 * x / (1 + x^2)) + 8 * cos(1.2*(k-1)) + sqrt(Q) * randn; %状态值
y = (x^2 / 20) + sqrt(R) * randn; %观测值
%扩展卡尔曼滤波器
%进行估计第一阶段的估计
e_x_estimate_1 = linear * e_x_estimate + 25 * e_x_estimate /(1+e_x_estimate^2) + 8 * cos(1.2*(k-1));
e_y_estimate = (e_x_estimate_1)^2/20; %这是根据k=1时估计值为1得到的观测值；只是这个由我估计得到的第24行的y也是观测值不过是由加了噪声的真实值得到的
%相关矩阵
e_A = linear + 25 * (1-e_x_estimate^2)/((1+e_x_estimate^2)^2);%传递矩阵
e_H = e_x_estimate_1/10; %观测矩阵
%估计的误差
e_p_estimate = e_A * e_P * e_A' + Q;
%扩展卡尔曼增益
e_K = e_p_estimate * e_H'/(e_H * e_p_estimate * e_H' + R);
%进行估计值的更新第二阶段
e_x_estimate_2 = e_x_estimate_1 + e_K * (y - e_y_estimate);
%更新后的估计值的误差
e_p_estimate_update = e_p_estimate - e_K * e_H * e_p_estimate;
%进入下一次迭代的参数变化
e_P = e_p_estimate_update;
e_x_estimate = e_x_estimate_2;
% 粒子滤波器
% 粒子滤波器
for i = 1 : N
p_xpartminus(i) = 0.5 * p_xpart(i) + 25 * p_xpart(i) / (1 + p_xpart(i)^2) + 8 * cos(1.2*(k-1)) + sqrt(Q) * randn; %这个式子比下面一行的效果好
% xpartminus(i) = 0.5 * xpart(i) + 25 * xpart(i) / (1 + xpart(i)^2) + 8 * cos(1.2*(k-1));
p_ypart = p_xpartminus(i)^2 / 20; %预测值
p_vhat = y - p_ypart;% 观测和预测的差
p_q(i) = (1 / sqrt(R) / sqrt(2*pi)) * exp(-p_vhat^2 / 2 / R); %各个粒子的权值
end
% 平均每一个估计的可能性
p_qsum = sum(p_q);
for i = 1 : N
p_q(i) = p_q(i) / p_qsum;%各个粒子进行权值归一化
end
% 重采样权重大的粒子多采点，权重小的粒子少采点, 相当于每一次都进行重采样；
for i = 1 : N
p_u = rand;
p_qtempsum = 0;
for j = 1 : N
p_qtempsum = p_qtempsum + p_q(j);
if p_qtempsum >= p_u
p_xpart(i) = p_xpartminus(j); %在这里xpart(i) 实现循环赋值；终于找到了这里！！！
break;
end
end
end
p_x_estimate = mean(p_xpart);
% p_x_estimate = 0;
% for i = 1 : N
% p_x_estimate =p_x_estimate + p_q(i)*p_xpart(i);
% end
%不敏卡尔曼滤波器
%采样点的选取存在x(i)
u_x_par = u_x_estimate;
for i = 2 : (u_symmetry_number+1)
u_x_par(i,:) = u_x_estimate + sqrt((u_symmetry_number+u_k) * u_P);
end
for i = (u_symmetry_number+2) : u_total_number
u_x_par(i,:) = u_x_estimate - sqrt((u_symmetry_number+u_k) * u_P);
end
%计算权值
u_w_1 = u_k/(u_symmetry_number+u_k);
u_w_N1 = 1/(2 * (u_symmetry_number+u_k));
%把这些粒子通过传递方程得到下一个状态
for i = 1: u_total_number
u_x_par(i) = 0.5 * u_x_par(i) + 25 * u_x_par(i)/(1+u_x_par(i)^2) + 8 * cos(1.2*(k-1));
end
%传递后的均值和方差
u_x_next = u_w_1 * u_x_par(1);
for i = 2 : u_total_number
u_x_next = u_x_next + u_w_N1 * u_x_par(i);
end
u_p_next = Q + u_w_1 * (u_x_par(1)-u_x_next) * (u_x_par(1)-u_x_next)';
for i = 2 : u_total_number
u_p_next = u_p_next + u_w_N1 * (u_x_par(i)-u_x_next) * (u_x_par(i)-u_x_next)';
end
% %对传递后的均值和方差进行采样产生粒子存在y(i)
% u_y_2obser(1) = u_x_next;
% for i = 2 : (u_symmetry_number+1)
% u_y_2obser(i,:) = u_x_next + sqrt((u_symmetry_number+k) * u_p_next);
% end
% for i = (u_symmetry_number + 2) : u_total_number
% u_y_2obser(i,:) = u_x_next - sqrt((u_symmetry_number+u_k) * u_p_next); % end
%另外存在y_2obser(i) 中；
for i = 1 :u_total_number
u_y_2obser(i,:) = u_x_par(i);
end
%通过观测方程得到一系列的粒子
for i = 1: u_total_number
u_y_2obser(i) = u_y_2obser(i)^2/20;
end
%通过观测方程后的均值y_obse
u_y_obse = u_w_1 * u_y_2obser(1);
for i = 2 : u_total_number
u_y_obse = u_y_obse + u_w_N1 * u_y_2obser(i);
end
%Pzz测量方差矩阵
u_pzz = R + u_w_1 * (u_y_2obser(1)-u_y_obse)*(u_y_2obser(1)-u_y_obse)';
for i = 2 : u_total_number
u_pzz = u_pzz + u_w_N1 * (u_y_2obser(i) - u_y_obse)*(u_y_2obser(i) - u_y_obse)';
end
%Pxz状态向量与测量值的协方差矩阵
u_pxz = u_w_1 * (u_x_par(1) - u_x_next)* (u_y_2obser(1)-u_y_obse)';
for i = 2 : u_total_number
u_pxz = u_pxz + u_w_N1 * (u_x_par(i) - u_x_next) * (u_y_2obser(i)- u_y_obse)';
end
%卡尔曼增益
u_K = u_pxz/u_pzz;
%估计量的更新
u_x_next_optimal = u_x_next + u_K * (y - u_y_obse);%第一步的估计值+ 修正值；
u_x_estimate = u_x_next_optimal;
%方差的更新
u_p_next_update = u_p_next - u_K * u_pzz * u_K';
u_P = u_p_next_update;
%进行画图程序
x_array = [x_array,x];
e_x_estimate_array = [e_x_estimate_array,e_x_estimate];
p_x_estimate_array = [p_x_estimate_array,p_x_estimate];
u_x_estimate_array = [u_x_estimate_array,u_x_estimate];
e_error(k,:) = abs(x_array(k)-e_x_estimate_array(k));
p_error(k,:) = abs(x_array(k)-p_x_estimate_array(k));
u_error(k,:) = abs(x_array(k)-u_x_estimate_array(k));
end
t = 0 : tf;
figure;
plot(t,x_array,'k.',t,e_x_estimate_array,'r-',t,p_x_estimate_array,'g--',t,u_x_estimate_array,'b:');
set(gca,'FontSize',10);
set(gcf,'color','White');
xlabel('时间步长');% lable --->label 我的神
ylabel('状态');
legend('真实值','EKF估计值','PF估计值','UKF估计值');
figure;
plot(t,x_array,'k.',t,p_x_estimate_array,'g--', t, p_x_estimate_array-1.96*sqrt(P), 'r:', t, p_x_estimate_array+1.96*sqrt(P), 'r:');
set(gca,'FontSize',10);
set(gcf,'color','White');
xlabel('时间步长');% lable --->label 我的神
ylabel('状态');
legend('真实值','PF估计值', '95% 置信区间');
%root mean square 平均值的平方根
e_xrms = sqrt((norm(x_array-e_x_estimate_array)^2)/tf);
disp(['EKF估计误差均方值=',num2str(e_xrms)]);
p_xrms = sqrt((norm(x_array-p_x_estimate_array)^2)/tf);
disp(['PF估计误差均方值=',num2str(p_xrms)]);
u_xrms = sqrt((norm(x_array-u_x_estimate_array)^2)/tf);
disp(['UKF估计误差均方值=',num2str(u_xrms)]);
% plot(t,e_error,'r-',t,p_error,'g--',t,u_error,'b:');
% legend('EKF估计值误差','PF估计值误差','UKF估计值误差');
t = 1 : tf;
figure;
plot(t,e_error,'r-',t,p_error,'g--',t,u_error,'b:');
set(gca,'FontSize',10);
set(gcf,'color','White');
xlabel('时间步长');% lable --->label 我的神
ylabel('状态');
legend('EKF估计值误差','PF估计值误差','UKF估计值误差');
% toc;
实验用例： 用模拟退火算法解决如下10个城市的TSP 问题，该问题最优解为691.2 opt f 。 表1 10个城市的坐标 城市 X 坐标 Y 坐标 城市 X 坐标 Y 坐标 3 0.4000 0.4439 8 0.8732 0.6536 编程实现 用MATLAB 实现模拟退火算法时，共编制了5个m 文件，分别如下 1、swap.m function [ newpath , position ] = swap( oldpath , number ) % 对 oldpath 进 行 互 换 操 作 % number 为 产 生 的 新 路 径 的 个 数 % position 为 对 应 newpath 互 换 的 位 置 m = length( oldpath ) ; % 城 市 的 个 数 newpath = zeros( number , m ) ; position = sort( randi( m , number , 2 ) , 2 ); % 随 机 产 生 交 换 的 位 置 for i = 1 : number newpath( i , : ) = oldpath ; % 交 换 路 径 中 选 中 的 城 市 newpath( i , position( i , 1 ) ) = oldpath( position( i , 2 ) ) ; newpath( i , position( i , 2 ) ) = oldpath( position( i , 1 ) ) ; end 2、pathfare.m function [ objval ] = pathfare( fare , path ) % 计 算 路 径 path 的 代 价 objval % path 为 1 到 n 的 排 列 ，代 表 城 市 的 访 问 顺 序 ； % fare 为 代 价 矩 阵 ， 且 为 方 阵 。 [ m , n ] = size( path ) ; objval = zeros( 1 , m ) ; for i = 1 : m for j = 2 : n objval( i ) = objval( i ) + fare( path( i , j - 1 ) , path( i , j ) ) ; end objval( i ) = objval( i ) + fare( path( i , n ) , path( i , 1 ) ) ; end
MCM战备历程3（模拟退火算法Matlab源程序）For glory 2007-02-03 11:20:04| 分类：数学建模 | 标签：学习|字号订阅 %模拟退火算法程序 T_max=input('please input the start temprature'); T_min=input('please input the end temprature'); iter_max=input('please input the most interp steps on the fit temp'); s_max=input('please input the most steady steps ont the fit temp'); T=T_max; load d:\address.txt; order1=randperm(size(address,1))';%生成初始解。 plot(address(order1,1),address(order1,2),'*r-') totaldis1=distance(address,order1); while T>=T_min iter_num=1; s_num=1; plot(T,totaldis1) hold on while iter_num
模拟退火算法模拟退火算法是一种通用的随机搜索算法，是局部搜索算法的扩展。它的思想是再1953 年由metropolis 提出来的，到1983 年由kirkpatrick 等人成功地应用在组合优化问题中。 模拟退火算法来源于固体退火原理，将固体加温至充分高，再让其徐徐冷却，加温时，固体内部粒子随温升变为无序状，内能增大，而徐徐冷却时粒子渐趋有序，在每个温度都达到平衡态，最后在常温时达到基态，内能减为最小。根据Metropolis 准则，粒子在温度T 时趋于平衡的概率为e- △ E/(kT)，其中E为温度T时的内能，AE为其改变量，k 为Boltzmann 常数。用固体退火模拟组合优化问题，将内能E模拟为目标函数值f ,温度T演化成控制参数t，即得到解组合优化问题的模拟退火算法：由初始解i和控制参数初值t开始，对当前解重复“产生新解-计算目标函数差T接受或舍弃”的迭代，并逐步衰减t值，算法终止时的当前解即为所得近似最优解，这是基于蒙特卡罗迭代求解法的一种启发式随机搜索过程。退火过程由冷却进度表(Cooli ng Schedule)控制，包括控制参数的初值t 及其衰减因子△ t、每个t值时的迭代次数L和停止条件S。 模拟退火算法新解的产生和接受可分为如下四个步骤：第一步是由一个产生函数从当前解产生一个位于解空间的新解；为便于后续的计算和接受，减少算法耗时，通常选择由当前新解经过简单地变换即可产生新解的方法，如对构成新解的全部或部分元素进行置换、互换等，注意到产生新解的变换方法决定了当前新解的邻域结构，因而对冷却进度表的选取有一定的影响。 第二步是计算与新解所对应的目标函数差。因为目标函数差仅由变换部分产生，所以目标函数差的计算最好按增量计算。事实表明，对大多数应用而言，这是计算目标函数差的最快方法。 第三步是判断新解是否被接受,判断的依据是一个接受准则，最常用的接受准则是Metropo1is准则：若厶t‘ <0 则接受S'作为新的当前解S,否则以概率exp(- △ t‘ /T) 接受S'作为新的当前解S。 第四步是当新解被确定接受时，用新解代替当前解，这只需将当前解中对应于产生新解时的变换部分予以实现，同时修正目标函数值即可。此时，当前解实现了一次迭代。 可在此基础上开始下一轮试验。而当新解被判定为舍弃时，
/* * 使用模拟退火算法(SA)求解TSP问题(以中国TSP问题为例) * 参考自《Matlab 智能算法30个案例分析》 * 模拟退火的原理这里略去，可以参考上书或者相关论文 * update: 16/12/11 * author:lyrichu * email:919987476@https://www.doczj.com/doc/e26558881.html, */ #include
第一章插值方法 1.1Lagrange插值 1.2逐步插值 1.3分段三次Hermite插值 1.4分段三次样条插值 第二章数值积分 2.1 Simpson公式 2.2 变步长梯形法 2.3 Romberg加速算法 2.4 三点Gauss公式 第三章常微分方程德差分方法 3.1 改进的Euler方法 3.2 四阶Runge-Kutta方法 3.3 二阶Adams预报校正系统 3.4 改进的四阶Adams预报校正系统 第四章方程求根 4.1 二分法 4.2 开方法 4.3 Newton下山法 4.4 快速弦截法 第五章线性方程组的迭代法 5.1 Jacobi迭代 5.2 Gauss-Seidel迭代 5.3 超松弛迭代 5.4 对称超松弛迭代 第六章线性方程组的直接法 6.1 追赶法 6.2 Cholesky方法 6.3 矩阵分解方法 6.4 Gauss列主元消去法
第一章插值方法 1.1Lagrange插值 计算Lagrange插值多项式在x=x0处的值. MATLAB文件：（文件名：Lagrange_eval.m）function [y0,N]= Lagrange_eval(X,Y,x0) %X,Y是已知插值点坐标 %x0是插值点 %y0是Lagrange插值多项式在x0处的值 %N是Lagrange插值函数的权系数 m=length(X); N=zeros(m,1); y0=0; for i=1:m N(i)=1; for j=1:m if j~=i; N(i)=N(i)*(x0-X(j))/(X(i)-X(j)); end end y0=y0+Y(i)*N(i); end 用法》X=[…];Y=[…]; 》x0= ; 》[y0,N]= Lagrange_eval(X,Y,x0) 1.2逐步插值 计算逐步插值多项式在x=x0处的值. MATLAB文件：（文件名：Neville_eval.m）function y0=Neville_eval(X,Y,x0) %X,Y是已知插值点坐标 %x0是插值点 %y0是Neville逐步插值多项式在x0处的值 m=length(X); P=zeros(m,1); P1=zeros(m,1); P=Y; for i=1:m P1=P; k=1; for j=i+1:m k=k+1;
模拟退火算法 模拟退火算法来源于固体退火原理，将固体加温至充分高，再让其徐徐冷却，加温时，固体内部粒子随温升变为无序状，内能增大，而徐徐冷却时粒子渐趋有序，在每个温度都达到平衡态，最后在常温时达到基态，内能减为最小。根据Metropolis准则，粒子在温度T时趋于平衡的概率为e-ΔE/(kT)，其中E为温度T时的内能，ΔE为其改变量，k为Boltzmann常数。用固体退火模拟组合优化问题，将内能E模拟为目标函数值f，温度T演化成控制参数t，即得到解组合优化问题的模拟退火算法：由初始解i和控制参数初值t开始，对当前解重复“产生新解→计算目标函数差→接受或舍弃”的迭代，并逐步衰减t值，算法终止时的当前解即为所得近似最优解，这是基于蒙特卡罗迭代求解法的一种启发式随机搜索过程。退火过程由冷却进度表(Cooling Schedule)控制，包括控制参数的初值t及其衰减因子Δt、每个t值时的迭代次数L和停止条件S。 模拟退火算法的模型 模拟退火算法可以分解为解空间、目标函数和初始解三部分。 模拟退火的基本思想: (1)初始化：初始温度T(充分大)，初始解状态S(是算法迭代的起 点)，每个T值的迭代次数L (2) 对k=1，……，L做第(3)至第6步： (3) 产生新解S′ (4) 计算增量Δt′=C(S′)-C(S)，其中C(S)为评价函数 (5) 若Δt′<0则接受S′作为新的当前解，否则以概率exp(-Δt′/T)
接受S′作为新的当前解. (6) 如果满足终止条件则输出当前解作为最优解，结束程序。终止条件通常取为连续若干个新解都没有被接受时终止算法。 (7) T逐渐减少，且T->0，然后转第2步。 算法对应动态演示图： 模拟退火算法新解的产生和接受可分为如下四个步骤： 第一步是由一个产生函数从当前解产生一个位于解空间的新解；为便于后续的计算和接受，减少算法耗时，通常选择由当前新解经过简单地变换即可产生新解的方法，如对构成新解的全部或部分元素进行置换、互换等，注意到产生新解的变换方法决定了当前新解的邻域结构，因而对冷却进度表的选取有一定的影响。 第二步是计算与新解所对应的目标函数差。因为目标函数差仅由变换部分产生，所以目标函数差的计算最好按增量计算。事实表明，对大多数应用而言，这是计算目标函数差的最快方法。 第三步是判断新解是否被接受,判断的依据是一个接受准则，最常用的接受准则是Metropo1is准则: 若Δt′<0则接受S′作为新的当前解S，否则以概率exp(-Δt′/T)接受S′作为新的当前解S。 第四步是当新解被确定接受时，用新解代替当前解，这只需将当前解中对应于产生新解时的变换部分予以实现，同时修正目标函数值即可。此时，当前解实现了一次迭代。可在此基础上开始下一轮试验。而当新解被判定为舍弃时，则
模拟退火算法及其Matlab 实现 模拟退火算法(Simulated Annealing algorithm ,简称SA )是柯克帕垂克(S. Kirkpatrick )于1982年受热力学中的固体退火过程与组合优化问题求解之间的某种“相似性”所启发而提出的,用于求解大规模组合优化问题的一种具有全局搜索 功能的随机性近似算法。与求解线性规划的单纯形法、Karmarkar 投影尺度法,求 解非线性规划的最速下降法、Newton 法、共轭梯度法,求解整数规划的分支定界法、割平面法等经典的优化算法相比,模拟退火算法在很大程度上不受制于优化问 题的具体形式和结构,具有很强的适应性和鲁棒性,因而也具有广泛的应用价值。 模拟退火算法源于对固体退火过程的模拟;采用Metropolis 接受准则;并用 一组称为冷却进度表的参数来控制算法进程,使得算法在多项式时间里给出一个近 似最优解。固体退火过程的物理现象和统计性质是模拟退火算法的物理背 景;Metropolis 接受准则使算法能够跳离局部最优的“陷阱”,是模拟退火算法能 够获得整体最优解的关键;而冷却进度表的合理选择是算法应用的关键。 1 物理退火过程 物理中的固体退火是先将固体加热至熔化,再徐徐冷却,使之凝固成规整晶体 的热力学过程。在加热固体时,固体粒子的热运动不断增加,随着温度的升高,粒子 与其平衡位置的偏离越来越大,当温度升至溶解温度后,固体的规则性被彻底破坏, 固体溶解为液体,粒子排列从较有序的结晶态转变为无序的液态,这个过程称为溶解。溶解过程的目的是消除系统中原先可能存在的非均匀状态,使随后进行的冷却 过程以某一平衡态为始点。溶解过程与系统的熵增过程相联系,系统能量也随温度 的升高而增大。 冷却时,液体粒子的热运动渐渐减弱,随着温度的徐徐降低,粒子运动渐趋有 序。当温度降至结晶温度后,粒子运动变为围绕晶体格点的微小振动,液体凝固成固体的晶态,这个过程称为退火。退火过程之所以必须“徐徐”进行,是为了使系统在每一温度下都达到平衡态,最终达到固体的基态(图1-1)。退火过程中系统的熵值
%%% 模拟退火算法源程序 % 此题以中国31省会城市的最短旅行路径为例： % clear;clc; function [MinD,BestPath]=MainAneal(pn) % CityPosition存储的为每个城市的二维坐标x和y； CityPosition=[1304 2312;3639 1315;4177 2244;3712 1399;3488 1535;3326 1556;3238 1229;... 4196 1044;4312 790;4386 570;3007 1970;2562 1756;2788 1491;2381 1676;... 1332 695;3715 1678;3918 2179;4061 2370;3780 2212;3676 2578;4029 2838;... 4263 2931;3429 1908;3507 2376;3394 2643;3439 3201;2935 3240;3140 3550;... 2545 2357;2778 2826;2370 2975]; figure(1); plot(CityPosition(:,1),CityPosition(:,2),'o') m=size(CityPosition,1);%城市的数目 % D = sqrt((CityPosition(:,ones(1,m)) - CityPosition(:,ones(1,m))').^2 + ... (CityPosition(:,2*ones(1,m)) - CityPosition(:,2*ones(1,m))').^2); path=zeros(pn,m); for i=1:pn path(i,:)=randperm(m); end iter_max=100;%i m_max=5;% Len1=zeros(1,pn);Len2=zeros(1,pn);path2=zeros(pn,m); t=zeros(1,pn); T=1e5; tau=1e-5; N=1; while T>=tau iter_num=1; m_num=1; while m_num 第8章 常微分方程初值问题数值解法 8.1 引言 8.2 欧拉方法 8.3 龙格-库塔方法 8.4 单步法的收敛性与稳定性 8.5 线性多步法 8.1 引 言 考虑一阶常微分方程的初值问题 00(,),[,],(). y f x y x a b y x y '=∈=（1.1） （1.2） 如果存在实数 ，使得 121212(,)(,).,R f x y f x y L y y y y -≤-?∈（1.3） 则称 关于 满足李普希茨(Lipschitz ）条件， 称为 的李普希茨常数（简称Lips.常数）. 0>L f y L f (参阅教材386页) 计算方法及MATLAB 实现 所谓数值解法，就是寻求解 在一系列离散节点 )(x y <<<<<+121n n x x x x 上的近似值 . ,,,,,121+n n y y y y 相邻两个节点的间距 称为步长. n n n x x h -=+1 如不特别说明，总是假定 为定数， ),2,1( ==i h h i 这时节点为 . ) ,2,1,0(0 =+=i nh x x n 初值问题(1.1),(1.2)的数值解法的基本特点是采取 “步进式”. 即求解过程顺着节点排列的次序一步一步地向前推进. 00(,),[,], (). y f x y x a b y x y '=∈= 描述这类算法，只要给出用已知信息 ,,,21--n n n y y y 计算 的递推公式. 1+n y 一类是计算 时只用到前一点的值 ，称为单步法. 1+n y n y 另一类是用到 前面 点的值 ， 1+n y k 11,,,+--k n n n y y y 称为 步法. k 其次，要研究公式的局部截断误差和阶，数值解 与 精确解 的误差估计及收敛性，还有递推公式的计算 稳定性等问题. n y )(n x y 首先对方程 离散化，建立求数值解的递推 公式. ),(y x f y =' function [f,T]=TSPSA(d,t0,tf) %TSP问题（货郎担问题，旅行商问题）的模拟退火算法通用malab源程序% f目标最优值,T最优路线,d距离矩阵,t0初始温度,tf结束温度 [m,n]=size(d); L=100*n; t=t0; pi0=1:n; min_f=0; for k=1：n-1 min_f=min_f+d(pi0(k),pi0(k+1)); end min_f=min_f+d(pi0(n),pi0(1)); p_min=pi0; while t>tf for k=1:L; kk=rand; [d_f,pi_1]=exchange_2(pi0,d); r_r=rand; if d_f<0 pi0=pi_1; elseif exp(d_f/t)>r_r pi0=pi_1; else pi0=pi0; end end f_temp=0; for k=1：n-1 f_temp=f_temp+d(pi0(k),pi0(k+1)); end f_temp=f_temp+d(pi0(n),pi0(1)); if min_f>f_temp min_f=f_temp; p_min=pi0; end t=0.87*t; end f=min_f; T=p_min; %aiwa要调用的子程序,用于产生新解 function [d_f,pi_r]=exchange_2(pi0,d) [m,n]=size(d); clear m; u=rand; u=u*(n-2); u=round(u); if u<2 u=2; end if u>n-2 u=n-2; end v=rand; v=v*(n-u+1); v=round(v); if v<1 v=1; end v=u+v; if v>n v=n; end pi_1(u)=pi0(v); pi_1(v)=pi0(u); if u>1 for k=1：u-1 pi_1(k)=pi0(k); end end if v>(u+1) for k=1：v-u-1 pi_1(u+k)=pi0(v-k); end end if v 用MATLAB实现结构可靠度计算 口徐华…朝泽刚‘u刘勇‘21 。 (【l】中国地质大学(武汉工程学院湖北?武汉430074; 12】河海大学土木工程学院江苏?南京210098 摘要:Matlab提供了各种矩阵的运算和操作,其中包含结构可靠度计算中常用的各种数值计算方法工具箱,本文从基本原理和相关算例分析两方面,阐述利用Matlab,编制了计算结构可靠度Matlab程.序,使得Matlab-语言在可靠度计算中得到应用。 关键词:结构可靠度Matlab软件最优化法 中图分类号:TP39文献标识码:A文章编号:1007-3973(200902-095-Ol 1结构可靠度的计算方法 当川概率描述结构的可靠性时,计算结构可靠度就是计算结构在规定时问内、规定条件F结构能够完成预定功能的概率。 从简单到复杂或精确稃度的不同,先后提出的可靠度计算方法有一次二阶矩方法、二次二阶矩方法、蒙特卡洛方法以及其他方法。一次■阶矩方法又分为。I-心点法和验算点法,其中验算点法足H前可靠度分析最常川的方法。 2最优化方法计算可靠度指标数学模型 由结构111n个任意分布的独立随机变量一,x:…以表示的结构极限状态方程为:Z=g(■.托…t=0,采用R-F将非正念变量当罱正态化,得到等效正态分布的均值o:和标准差虹及可靠度指标B,由可靠度指标B的几何意义知。o;辟 开始时验算点未知,把6看成极限状态曲面上点P(■,爿:---37,的函数,通过优化求解,找到B最小值。求解可靠皮指标aJ以归结为以下约束优化模型: rain睁喜t华,2 s.,.Z=g(工i,x2’,…,工:=0 如极限状态方栉巾某个变最(X。可用其他变量表示,则上述模型jfIJ‘转化为无约束优化模型: 。。B!:手f生丛r+阻:坚:坠:盐尘}二剐 t∞oY?’【叫,J 3用MATLAB实现结构可靠度计算 3.1Matlab简介 Matlab是++种功能强、效率高、便.丁.进行科学和工程计算的交互式软件包,汇集了人量数学、统计、科学和工程所需的函数,MATI.AB具有编程简甲直观、用户界mf友善、开放性强等特点。将MATLAB用于蒙特卡罗法的一个显著优点是它拥有功能强大的随机数发生器指令。 3.2算例 3.2.I例:已知非线形极限状态方程z=g(t r'H=567f r-0.5H2=0’f、r服从正态分布。IIf=0.6,o r=0.0786;la|_ 2.18,o r_0.0654;H服从对数正态分布。u H= 3218,O。 =0.984。f、r、H相互独立,求可靠度指标B及验算点(,,r’,H‘。 解:先将H当量正念化:h=ln H服从正态分布,且 ,‘-““了:等专虿’=,。49?口二-、『五ir面_。。3 模拟退火算法 基于模拟退火算法的TSP问题求解具体步骤如下： 1）随机产生一个初始解path（作为当前最优路径），计算目标函数值pathfare(fare,path)=e0，并设置初始温度t0，内循环终止方差istd，外循环终止方差ostd降温系数lam，温度更新函数tk=lam*tk-1，并令k=1，输入各城市坐标coord，计算城市间的距离fare。 2）根据控制参数更新函数来控制温度的下降过程，给定最大循环步数iLK，设置循环计数器的初始值in=1。 3）对当前的最优路径作随机变动，产生一个新路径newpath，计算新路径的目标函数值pathfare(fare,newpath)=e1和目标函数值的增量e1-e04）根据Metropolis准则，如果增量（e1-e0）<0，则接受新产生的路径newpath作为当前最优路径；如果（e1-e0）>=0，则以公式（1）来决定新路径newpath是否代替path。rand()随机产生一个在[0,1]之间的随机数。exp[-(e1-e0)/t]>rand() 4）如果目标函数值小于istd，则直接跳出内循环。 5）如果in 1. distance.m function [fare]=distance(coord) %coord为各城市的坐标 %fare为城市间的距离矩阵 [~,m]=size(coord);%m为城市的个数 fare=zeros(m); for i=1:m%外层为行 for j=1:m%内层为列 fare(i,j)=(sum((coord(:,i)-coord(:,j)).^2))^0.5; fare(j,i)=fare(i,j);%距离矩阵对称 end end 2. myplot.m function []=myplot(path,coord,pathfar) %做出路径的图形 %path为要做图的路径,coord为各个城市的坐标 %pathfar为路径path对应的费用 len=length(path); clf; hold on; 作者：夏云木子 1、 >> syms re(x) re(y) re(z) >> input('计算相对误差：'),re(x)=10/1991,re(y)=0.0001/1.991,re(y)=0.0000001/0.0001991 所以可知re(y)最小，即y精度最高 2、 >> format short,A=sqrt(2) >> format short e,B=sqrt(2) >> format short g,C=sqrt(2) >> format long,D=sqrt(2) >> format long e,E=sqrt(2) >> format long g,F=sqrt(2) >> format bank,H=sqrt(2) >> format hex,I=sqrt(2) >> format +,J=sqrt(2) >> format,K=sqrt(2) 3、 >> syms A >> A=[sqrt(3) exp(7);sin(5) log(4)];vpa(pi*A,6) 4、1/6251-1/6252=1/6251*6252 5、（1）1/（1+3x）-(1-x)/(1+x)=x*(3*x-1)/[(1+3*x)*(1+x)] (2) sqrt(x+1/x)-sqrt(x-1/x)=2/x/[sqrt(x-1/x)+sqrt(x+1/x)] (3) log10(x1)-log(x2)=log10(x1/x2) (4) [1-cos(2*x)]/x =x^2/factorial(2)-x^4/factorial(4)+x^6/factorial(6)-… 基于matlab的模拟退火法 编写一个matlab的程序用模拟退火法求函数最优解 function [xo,fo] = Opt_Simu(f,x0,l,u,kmax,q,TolFun) % 模拟退火算法求函数f(x)的最小值点，且l <= x <= u % f为待求函数，x0为初值点，l，u分别为搜索区间的上下限，kmax为最大迭代次数 % q为退火因子，TolFun为函数容许误差 %%%%算法第一步根据输入变量数，将某些量设为缺省值 if nargin < 7 TolFun = 1e-8; end if nargin < 6 q = 1; end if nargin < 5 kmax = 100; end %%%%算法第二步，求解一些基本变量 N = length(x0); %自变量维数 x = x0; fx = feval(f,x); %函数在初始点x0处的函数值 xo = x; fo = fx; %%%%%算法第三步，进行迭代计算，找出近似全局最小点 for k =0:kmax Ti = (k/kmax)^q; mu = 10^(Ti*100); % 计算mu dx = Mu_Inv(2*rand(size(x))-1,mu).*(u - l);%步长dx x1 = x + dx; %下一个估计点 x1 = (x1 < l).*l +(l <= x1).*(x1 <= u).*x1 +(u < x1).*u; %将x1限定在区间[l,u]上 fx1 = feval(f,x1); df = fx1- fx; if df < 0||rand < exp(-Ti*df/(abs(fx) + eps)/TolFun) %如果fx1 1、Newdon迭代法求解非线性方程 function [x k t]=NewdonToEquation(f,df,x0,eps) %牛顿迭代法解线性方程 %[x k t]=NewdonToEquation(f,df,x0,eps) %x:近似解 %k:迭代次数 %t:运算时间 %f:原函数，定义为内联函数 ?:函数的倒数，定义为内联函数 %x0:初始值 %eps:误差限 % %应用举例： %f=inline('x^3+4*x^2-10'); ?=inline('3*x^2+8*x'); %x=NewdonToEquation(f,df,1,0.5e-6) %[x k]=NewdonToEquation(f,df,1,0.5e-6) %[x k t]=NewdonToEquation(f,df,1,0.5e-6) %函数的最后一个参数也可以不写。默认情况下，eps=0.5e-6 %[x k t]=NewdonToEquation(f,df,1) if nargin==3 eps="0".5e-6; end tic; k=0; while 1 x="x0-f"(x0)./df(x0); k="k"+1; if abs(x-x0) < eps || k >30 break; end x0=x; end t=toc; if k >= 30 disp('迭代次数太多。'); x="0"; t="0"; end 2、Newdon迭代法求解非线性方程组 function y="NewdonF"(x) %牛顿迭代法解非线性方程组的测试函数 %定义是必须定义为列向量 y(1,1)=x(1).^2-10*x(1)+x(2).^2+8; y(2,1)=x(1).*x(2).^2+x(1)-10*x(2)+8; return; function y="NewdonDF"(x) %牛顿迭代法解非线性方程组的测试函数的导数 y(1,1)=2*x(1)-10; y(1,2)=2*x(2); y(2,1)=x(2).^+1; y(2,2)=2*x(1).*x(2)-10; return; 以上两个函数仅供下面程序的测试 function [x k t]=NewdonToEquations(f,df,x0,eps) %牛顿迭代法解非线性方程组 %[x k t]=NewdonToEquations(f,df,x0,eps) %x:近似解 %k:迭代次数 %t:运算时间 %f:方程组（事先定义） ?:方程组的导数（事先定义） %x0:初始值 %eps:误差限 % %说明：由于虚参f和df的类型都是函数，使用前需要事先在当前目录下采用函数M文件定义% 另外在使用此函数求解非线性方程组时，需要在函数名前加符号“@”，如下所示 % %应用举例： %x0=[0,0];eps=0.5e-6; %x=NewdonToEquations(@NewdonF,@NewdonDF,x0,eps) %[x k]=NewdonToEquations(@NewdonF,@NewdonDF,x0,eps) %[x k t]=NewdonToEquations(@NewdonF,@NewdonDF,x0,eps) %函数的最后一个参数也可以不写。默认情况下，eps=0.5e-6 %[x k t]=NewdonToEquations(@NewdonF,@NewdonDF,x0,eps) % maxpop给定群体规模% pop群体 % newpop种群 %t0初始温度 function [codmin,finmin]=fc0(cc,v0,t0) N=length(cc(1,:)); %定群体规模 if N>50 maxpop=2*N-20; end if N<=40 maxpop=2*N; end %产生初始群体 pop=zeros(maxpop,N); pop(:,1)=v0; finmin=inf; codmin=0; for i=1:maxpop Ra=randperm(N); Ra(find(Ra==v0))=Ra(1); pop(i,:)=Ra; end t=t0; while t>0 %用模拟退火产生新的群体pop=fc1(maxpop,pop,N,cc,v0,t); %转轮赌选择种群 f=zeros(1,maxpop); for i=1:maxpop for j=1:N-1 x=pop(i,j); y=pop(i,j+1); fo1=cc(pop(i,j),pop(i,j+1)); f(i)=f(i)+fo1; end f(i)=f(i)+cc(pop(i,1),pop(i,N)); end fmin=min(f); for i=1:maxpop if fmin==inf&f(i)==inf end if fmin~=inf|f(i)~=inf dd=fmin-f(i); end ftk(i)=exp(dd/t); end [fin1,cod]=sort(-ftk); fin=abs(fin1); %f(cod(1)) if f(cod(1)) 信息系统分析与设计作业 层次分析法确定绩效评价权重在matlab中的实现 小组成员：孙高茹、王靖、李春梅、郭荣1 程序简要概述 编写程序一步实现评价指标特征值lam、特征向量w以及一致性比率CR的求解。 具体的操作步骤是：首先构造评价指标，用专家评定法对指标两两打分，构建比较矩阵，继而运用编写程序实现层次分析法在MATLAB中的应用。 通过编写MATLAB程序一步实现问题求解，可以简化权重计算方法与步骤，减少工作量，从而提高人力资源管理中绩效考核的科学化电算化。 2 程序在matlab中实现的具体步骤 function [w,lam,CR] = ccfx(A) %A为成对比较矩阵，返回值w为近似特征向量 % lam为近似最大特征值λmax，CR为一致性比率 n=length(A(:,1)); a=sum(A); B=A %用B代替A做计算 for j=1:n %将A的列向量归一化 B(:,j)=B(:,j)./a(j); end s=B(:,1); for j=2:n s=s+B(:,j); end c=sum(s);%计算近似最大特征值λmax w=s./c; d=A*w lam=1/n*sum((d./w)); CI=(lam-n)/(n-1);%一致性指标 RI=[0,0,0.58,0.90,1.12,1.24,1.32,1.41,1.45,1.49,1.51];%RI为随机一致 性指标 CR=CI/RI(n);%求一致性比率 if CR>0.1 disp('没有通过一致性检验'); else disp('通过一致性检验'); end end 3 案例应用 我们拟构建公司员工绩效评价分析权重，完整操作步骤如下： 3.1构建的评价指标体系 我们将影响员工绩效评定的指标因素分为：打卡、业绩、创新、态度与品德。 3.2专家打分，构建两两比较矩阵 A = 1.0000 0.5000 3.0000 4.0000 2.0000 1.0000 5.0000 3.0000 0.3333 0.2000 1.0000 2.0000 0.2500 0.3333 0.5000 1.0000 3.3在MATLAB中运用编写好的程序实现 直接在MATLAB命令窗口中输入 [w,lam,CR]=ccfx(A) 继而直接得出 d = 1.3035 2.0000 0.5145 0.3926 w = 0.3102 0.4691 0.1242 0.0966 lam =4.1687 计算方法作业（作者：夏云木子） 1、help linspace type linspace 2、a1=[5 12 47;13 41 2;9 6 71];a2=[12 9;6 15;7 21];B=a1*a2, C=a1(:,1:2).*a2, D=a1.^2, E=a1(:).^2 3、a1=[5 12 47;13 41 2;9 6 71];a2=[12 9;6 15;7 21];a1(4:5,1:3)=a2.';a1([4 5],:)=a1([5 4],:);b1=a1 c1=b1(4,1),c2=b1(5,3),D=b1(3:4,:)*a2 4、a1=[5 12 47;13 41 2;9 6 71]; E=eye(3,3); S = a1 + 5*a1' - E, S1=a1^3-rot90(a1)^2+6*E 5、a1=[5 12 47;13 41 2;9 6 71];s=5;A=s-a1,B=s*a1,C=s.*a1,D=s./a1,E=a1./s 6、c=[1 2 3 4;5 6 7 8;9 10 11 12;13 14 15 16];A=c^-4,B=(c^3)^-1,C=(3*c+5*c^-1)/5 7、a=[1 i 3;9i 2-i 8;7 4 8+i];A=a.' 8、abc=[-2.57 8.87;-0.57 3.2-5.5i];m1=sign(abc),m2=round(abc),m3=floor(abc) Sign为符号函数，round表示四舍五入取整，floor表示舍去小数部分取整 9、x=[1 4 3 2 0 8 10 5]';y=[8 0 0 4 2 1 9 11]';A=dot(x,y) 10、a=[3.82 5.71 9.62];b=[7.31 6.42 2.48];A=dot(a,b),B=cross(a,b) 11、P=[5 7 8 0 1];Pf=poly(P);Px=poly2str(Pf,'x') 12、P=[3 0 9 60 0 -90];K1=polyval(P,45),K2=polyval(P,-123),K3=polyval(P,579) 13、P1=[13 55 0 -17 9];P2=[63 0 26 -85 0 105];PP=conv(P1,P2);P1P2=poly2str(PP,'x'),[Q,r]=deconv(P2,P1)0计算方法及MATLAB实现简明讲义课件PPS8-1欧拉龙格法
模拟退火算法求解TSP问题Matlab源码
用MATLAB实现结构可靠度计算.
M模拟退火最短距离问题 Matlab代码 亲测完美运行
计算方法及其MATLAB实现第二章作业
基于matlab的模拟退火法
matlab用于计算方法的源程序
遗传模拟退火算法matlab通用源程序
层次分析法计算权重在matlab中的实现
计算方法及其MATLAB实现第一章作业
相关主题
文本预览