变分法
变分法笔记(1)——古典变分问题的例子
上面作者的专栏
哈密顿力学专栏
变分法介绍专栏
matlab画带约束的三维图
搜了半天也没找到能 matlab画带约束的三维图
用约束constraint。照着这个修改吧
x=-2:0.1:10;
y=-2:0.1:10;
[X1,X2]=meshgrid(x,y);
constraint1 = -X1<=0;
constraint2 = -X2<=0;
constraint3 = X1+X2-11<=0;
constraint4 = X1-6<=0;
constraint5 = X2-8<=0;
X1(~constraint1) = NaN;
X2(~constraint1) = NaN;
X1(~constraint2) = NaN;
X2(~constraint2) = NaN;
X1(~constraint3) = NaN;
X2(~constraint3) = NaN;
X1(~constraint4) = NaN;
X2(~constraint4) = NaN;
X1(~constraint5) = NaN;
X2(~constraint5) = NaN;
f = X1.^2+X2.^2-X1.*X2-10*X1-4*X2+60;
surf(X1,X2,f);
最速下降法
syms x1 x2 s; %声明符号变量
% f1 = 8*x1^2+4*x1*x2+5*x2^2;
%f1 = x1^4-2*x1^2*x2-2*x1*x2+x1^2+2*x2^2+4.5*x1-4*x2+4;
f1 = x1.^2+x2.^2-x1.*x2-10*x1-4*x2+60;
%f1 = 100 * (x2 - x1).^2 + (1-x1).^2;
k=descent(f1,x1,x2,s,[2,2],10^(-3));
result_string=sprintf('在 %d 次迭代后求出极小值\n',k);%在迭代多少次之后求出极小值
disp(result_string);
function k = descent(f,x1,x2,s,x0,e)
k = 0;
grad_f = [diff(f,x1) diff(f,x2)];
delta = subs(grad_f,[x1,x2],[x0(1),x0(2)]);
step=1;
i = x0; %当前点i
while norm(delta) > e
k = k + 1;
x_next = [i(1),i(2)] - s* delta/norm(delta);
f_val = subs(f,[x1,x2],[x_next(1),x_next(2)]);
step = abs(double(solve(diff(f_val,s))));
step = step(1);
%x(k+1)点
i = double([i(1),i(2)] - step * delta/norm(delta));
%x(k+1)点梯度值
delta = subs(grad_f,[x1,x2],[i(1),i(2)]);
%函数值
f_value = double(subs(f,[x1,x2],[i(1),i(2)]));
%输出迭代计算过程
result_string=sprintf('k=%d, x1=%.6f, x2=%.6f, step=%.6f f(x1,x2)=%.6f',...
k,i(1),i(2),step,f_value);
disp(result_string);
end
end
function [x,val,k] = grad(fun,gfun,x0)
%功能：用最速下降法求解无约束问题 minif（x）
%输入：fun，gfun分别是目标函数和梯度，x0是初始点
%输出：x,val分别是近似最优值和最优值，k是迭代次数
maxk=5000; %最大迭代次数
rho=0.5;
sigma=0.4;
k=0;
e=0.005; %精度
while(k<maxk)
g=feval(gfun,x0); %计算梯度
d=-g;
if(norm(d)<e),break;end
%用Amrijo搜索技术确定步长
m=0;mk=0;
while(m<20) %最大迭代次数
if(feval(fun,x0+rho^m*d)<feval(fun,x0)+sigma*rho^m*g'*d)
mk=m;
break;
else
m=m+1;
end
end
x0=x0+d*rho^mk;
k=k+1;
end
x=x0;
val=feval(fun,x0);
end
clear all
clc
%利用grad函数求解 minif(x)=100*(x1^2-x2)^2+(x1-1)^2
x0=[10 10]';
[x,val,k]=grad('fun','gfun',x0);
disp(['最优解：x = '])
disp(x)
disp(['此时: f(x) = ',num2str(val)])
function g=gfun(x)
%目标函数的梯度
%g = [8*x(1)-2*x(1)*x(2), 2*x(2)-x(1)^2]'
%g=[400*x(1)*(x(1)^2-x(2))+2*(x(1)-1),-200*(x(1)^2-x(2))]';
g = [4*x(1)^3-4*x(1)*x(2)-2*x(2)+2*x(1)+4.5,-2*x(1)^2-2*x(1)+4*x(2)-4]';
% g = [16*x(1)+4*x(2),4*x(1)+10*x(2)]';
end
function f= fun(x)
%目标函数
%f=4*x(1)^2+x(2)^2-x(1)^2*x(2);
%f=100*(x(1)^2-x(2))^2+(x(1)-1)^2;
f = x(1)^4-2*x(1)^2*x(2)-2*x(1)*x(2)+x(1)^2+2*x(2)^2+4.5*x(1)-4*x(2)+4;
% f = 8*x(1)^2+4*x(1)*x(2)+5*x(2)^2;
end
单纯形法
单纯形优化法
单纯形的优化设计
matlab fminsearch 函数
单纯形代码 matlab
单纯形探究与改进
基于单纯形的最优化方法
作业：改进单纯形法并验证效果（无约束问题）
syms x1 x2
y = 4*(x1-5)^2+(x2-6)^2
a = [8,11]; % f = 65
b = [8,9]; % f = 45
c = [10 11]; % f = 125 最差点
% a = [4.706 8.353]; % f = 5.8824
% b = [8 9];
% c = [8 11]; % 最差点
% a = [4.9423 6.1786]; % f = 0.0452
% b = [4.706 8.353];
% c = [8 9]; % 最差点
% a = [4.9423 6.1786];
% c = [4.9619 5.9474];
% b = [5.0321 6.0712];
% 判断三个函数值 b<a<c
r1 = [(a(1)+b(1))/2 (a(2)+b(2))/2] % 形心点
r2 = c % 最差点
p = polyfit(r1,r2,1)
x3 = p(1)*x1+p(2)
f = subs(y,x2,x3) % 将直线代入 y
xishu = coeffs(f) % 得到方程系数
minx = double(-xishu(2)/(2*xishu(3))) %-b/2a 求函数最小值对应的x
minx2 = double(subs(x3,x1,minx))
fmin = double(subs(f,x1,minx)) % 函数最小值
自己写的代码，没有设置循环，每次运行得出的最优点代替原来的最差点，再修改 abc 矩阵。只针对直线带进去后为一元二次函数，所以用-2a 分之 b 求最小值。如果是高次的还得另外写求最小值。
书上迭代四次到 0.25
我的代码能到 0.04，但是0.04 再往下叠最小值就越来越大了，可能需要写出循环代码才能看什么问题。
还算了个课后题，两次就能出结果。
有约束问题
f = x(1).^2+x(2).^2-x(1).*x(2)-10*x(1)-4*x(2)+60;
用自己改的方法求解第五章课后的带约束题目：
最速下降法 无约束时 起始点[2 2] 精度 0.005：
fmincon函数无约束时
有约束时，fmincon 函数
function f = fminx(x)
f = x(1).^2+x(2).^2-x(1).*x(2)-10*x(1)-4*x(2)+60;
end
a = [1 0; -1 0; 0 1; 0 -1; 1 1]
b = [6;0;8;0;11]
[x,fval,exitflag] = fmincon(@(x) fminx(x),[0,0],a,b)
% 对于fmincon函数，其exitflag参数中的数字：
% 1、一阶最优性条件满足容许范围
% 2、X的变化小于容许范围
% 3、目标函数的变化小于容许范围
% 4、重要搜索方向小于规定的容许范围并且约束违背小于options.TolCon
% 5、重要方向导数小于规定的容许范围并且约束违背小于options.TolCon
% 0、到达最大迭代次数或到达函数评价
% -1、算法由输出函数终止
% -2、无可行点
这里用最速下降法的代码修改了一下，将 迭代获得的x0修改为该点在直线上的垂足，因为得到的结果很难满足 x1 和 x2 的上下限，所以就先迭代 50 次的结果，保存在 x3，x4上，对应的函数值保存在 result 中，但是这三个值不能保存在工作区，只能在命令窗口出现，于是使用 save，保存为文件然后再打开进行后一步操作。
function [x,val,k] = grad(fun,gfun,x0)
%fun 和 gfun 函数在前面
%输入：fun，gfun分别是目标函数和梯度，x0是初始点
%输出：x,val分别是近似最优值和最优值，k是迭代次数
maxk=50; %最大迭代次数
rho=0.5;
sigma=0.4;
k=0;
e=0.005; %精度
while(k<maxk)
g=feval(gfun,x0); %计算梯度
d=-g;
if(norm(d)<e)&(x0(1)<=6)&(x0(2)<=8),break;end
%用Amrijo搜索技术确定步长
m=0;mk=0;
while(m<20) %最大迭代次数
if(feval(fun,x0+rho^m*d)<feval(fun,x0)+sigma*rho^m*g'*d)
mk=m;
break;
else
m=m+1;
end
% if (x0(1)<=6)&(x0(2)<=8)
% else
% break;
% end
end
x0=x0+d*rho^mk;
x2 = [(11+x0(1)-x0(2))/2 11-( 11+x0(1)-x0(2) )/2 ];
x0 = x2
k=k+1
x3(k) = x0(1);
x4(k) = x0(2);
result(k)=feval(fun,x0);
end
x=x0;
val=feval(fun,x0);
result1 = result
save result1
end
得到 result 后，先获取满足 x3x4 的值，获得逻辑值矩阵，再与结果点乘，就能剩下筛选过的函数值，然后获取次小值（因为矩阵最小值是 0，所以要获取最小函数值得除了 0 以外的最小），然后再获取函数值对应的 x 值。
lo = (x3<=6)&(x4<=8)
valu = result.*lo
a=min(valu(find(valu-min(valu))))
x3(find(valu==a))
x4(find(valu==a))
zoutendijk法求解线性约束的非线性目标函数的最优解
function [X0,f_val]=zoutendijk(A,b,x0,Aeq,beq)
%自定义函数diff_val(x0)作用是求所给函数在x0出的偏导数
%自定义函数fval(x0)作用是求所给函数在x0出的函数值
function [X0,f_val]=zoutendijk(A,b,x0,Aeq,beq)
%自定义函数diff_val(x0)作用是求所给函数在x0出的偏导数
%自定义函数fval(x0)作用是求所给函数在x0出的函数值
format long;
eps=1.0e-6;
x0=transpose(x0); %刚开始给的x0为行向量
[f,x]=func;
sz=length(x0);
[m,n]=size(A);
%把约束方程系数矩阵A分解为A1,A2，其中A1为起作用约束
for k=1:1:100
A1=A;
A2=A;
b1=b;
b2=b;
for i=m:-1:1
if abs(A2(i,:)*x0-b2(i,:)) < 0.1
A2(i,:)=[];
b2(i,:)=[];
end
end
for i=m:-1:1
if abs(A1(i,:)*x0-b1(i,:))>=0.1
A1(i,:)=[];
b1(i,:)=[];
end
end
A1;
A2;
b1;
b2;
i2=rank(A2);
AE=[A1;Aeq];
[i1,j1]=size(AE);
r=rank(AE);
if r<i1
disp('行不满秩')
return
end
if i2==0
disp('无效')
return
end
%求解线性规划问题得到可行下降方向d0
s=diff_val(x0');
c=double(s);
lb=-1*ones(sz,1);
ub=ones(sz,1);
k1=length(b1);
k2=length(beq);
p=zeros(k1,1);
q=zeros(k2,1);
% [d0,mn,m1,m2,m3]=linprog(c,A1,p,Aeq,q,lb,ub);
[d0,mn]=linprog(c,A1,p,Aeq,q,lb,ub);
%[x,fval]= linprog(f,A,b,Aeq,beq,lb,ub)
% x返回的是决策向量的取值，fval返回的是目标函数最优解
% A，b：线性不等式约束，Aeq，beq：线性等式约束，ub，lb：决策向量上下界向量
d0;mn;
df=abs(s*d0);
%该点处的可行下降方向与该点梯度乘积判定终止迭代
if df<0.1
disp( '最优解为:')
x0
f_val=fval(x0)
k
return
else
%进行一维搜索，求f(x(k+1))的最小值
b_=b2-A2*x0;
d_=A2*d0;
[dh,dl]=size(d_);
ul=1;
for i=1:1:dh
if d_(i,:)>=0
u=1;
else
u=0;
end
ul=ul*u;
end
ul;b_;d_;
vmax=inf;
if ul==0
vmax=inf;
else
for i=1:1:dh
if d_(i,:)>0
v=b_(i,:)/d_(i,:);
if v<vmax
vmax=v;
end
end
end
end
end
vmax;
h=fmin(x0,d0,vmax);
a=x0+h*d0;
f_val=fval(a);
x0=x0+h*d0; %新的可行点
% '****************'
X0=x0
f_val=fval(x0)
k
% '****************'
end
end
% A = [1 2 1; -1 0 0; 0 -1 0; 0 0 -1];
% b = [4 0 0 0]';
% x0 = [0 0 0];
% Aeq = [];
% beq = [];
% label = 1;
% zoutendijk(A,b,x0,Aeq,beq)
% 输入
% https://blog.csdn.net/weixin_43290709/article/details/122441446
A = [1 2 1; -1 0 0; 0 -1 0; 0 0 -1];
b = [4 0 0 0]';
x0 = [0 0 0];
Aeq = [];
beq = [];
zoutendijk(A,b,x0,Aeq,beq)
% h=fmin(x0,d0,vmax)
% d0 = [1;1;1]
% x0 = [0;0;0]
% vmax = 1
% A = [2 -1; 1 1; -1 0; 0 -1];
% b = [1 2 0 0]';
% x0 = [0 0];
% Aeq = [];
% beq = [];
% zoutendijk(A,b,x0,Aeq,beq)
function [f,x]=func
%设置目标函数
syms x1 x2 x3;
f=x1^2+2*x2^2+3*x3^2+x1*x2-2*x1*x3+x2*x3-4*x1-6*x2;
%f = (x1-1)^2+(x2-2)^2+1;
x=[x1,x2,x3];
% x=[x1,x2];
end
function f_val=fval(x0)
%求目标函数值
x0=transpose(x0);
[f,x]=func;
f_val=double(subs(f,x,x0));
end
function h=fmin(x0,d0,vmax)
%求函数最小值
[f,x]=func;
syms h;
a=x0+h*d0;
f_val=inline(subs(f,x,a')); % a为 3*1 x 需要 转置
if vmax==inf
min_h=fminbnd(f_val,0,10000);
else
min_h=fminbnd(f_val,0,vmax);% x=fminbnd(fun,x1,x2)
% 返回目标函数f(x)在区间（x1,x2）上的极小值
end
h=min_h;
end
5914
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
