You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
mycode/NLP.m

136 lines
6.5 KiB

This file contains ambiguous Unicode characters!

This file contains ambiguous Unicode characters that may be confused with others in your current locale. If your use case is intentional and legitimate, you can safely ignore this warning. Use the Escape button to highlight these characters.

%非线性规划
%P42例题
%求解
clc,clear,format long g,close all;
syms x1 x2;%符号变量(用于定义函数、求解方程、计算导数和积分)
f=(339-0.01*x1-0.003*x2)*x1+(399-0.004*x1-0.01*x2)*x2-(195*x1+225*x2+400000);
f=simplify(f);
f1=diff(f,x1),f2=diff(f,x2);
[x10,x20]=solve(f1,f2);%求驻点即使f1,f2为0的点
x10=round(double(x10));%四舍五入
x20=round(double(x20));
f0=subs(f,{x1,x2},{x10,x20});
f0=double(f0);
subplot(121);
fmesh(f,[0,10000,0,10000]);%3维
%xlabel('x1','Interpreter','latex');
%ylabel('y1','Interpreter','latex');
xlabel('x1');
ylabel('y1');
subplot(122);
L=fcontour(f,[0,10000,0,10000]);
contour(L.XData,L.YData,L.ZData,'ShowText','on');
xlabel('x1');
ylabel('y1');
p1=339-0.01*x10-0.003*x20;
p2=399-0.004*x10-0.01*x20;
c=195*x10+225*x20+400000;
rate=f0/c;
%灵敏性检验 分析 最优解 对 19村彩电价格下降幅度 的灵敏性
clc,clear,close all,format long g;
syms x1 x2 a;
X1=4735,X2=7043,AA=0.01,FF=553641;
f=(339-a*x1-0.003*x2)*x1+(399-0.004*x1-0.01*x2)*x2-(195*x1+225*x2+400000);
f1=diff(f,x1),f2=diff(f,x2);
[x10,x20]=solve(f1,f2);
subplot(121);
fplot(x10,[0.002,0.02]);
xlabel('a');
ylabel('x1');
subplot(122);
fplot(x20,[0.002,0.02]);
xlabel('a');
ylabel('x2');
dx1=diff(x10,a);
dx10=subs(dx1,a,0.01),dx10=double(dx10);
sx1a=dx10*AA/X1;
dx2=diff(x20,a);
dx20=subs(dx2,a,0.01),dx20=double(dx20);
sx2a=dx20*AA/X2;
F=subs(f,[x1,x2],[x10,x20]);
F=simplify(F);
figure,fplot(F,[0.002,0.02]);
xlabel('a');
ylabel('profit');
dpro=diff(F,a);
dpro0=subs(F,a,0.01);
sya=dpro0*AA/FF;
%当a提高10%
f3=subs(f,[x1,x2,a],[X1,X2,0.011]),f3=double(f3);
f4=subs(F,[a],0.011);
%相对误差
delta=(f4-f3)/f4;
delta=double(delta);
%相对误差很小说明在19寸的价格下降幅度在一定范围内发生变化时仍然按照不变时的最优解来确定生产方案仍然可以获得几乎最优的利润。
%使用fmincon的简单例子
%蒙特卡罗确定初始值
clc,clear,close all;
n=10^7;
x1=unifrnd(0,12,[1,n]);
x2=unifrnd(0,12,[1,n]);
fmin=+inf;
for i=1:n
x=[x1(i),x2(i)];
if ((x1(i)-1)^2-x2(i)-6<=0) & (-2*x1(i)+3*x2(i)-6<=0)
temp=x(1)^2+x(2)^2-x(1)*x(2)-2*x(1)-5*x(2);
if temp<fmin
fmin=temp;
x0=x;
end
end
end
%[x,fval]= fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon,options)
A=[-2,3];
b=6;
[x,fval]=fmincon(@fun1,x0,A,b,[],[],[],[],@nonlfun1);
%option=optimoptions('fmincon','Algorithm','sqp');
%[x,fval]=fmincon(@fun1,x0,A,b,[],[],[],[],@nonlfun1,option);
disp('Thw optimal solution is:');
disp(x);
disp('The minimum value is:');
disp(fval);
figure;
fmesh(@(x1, x2) arrayfun(@(x1, x2) fun1([x1, x2]), x1, x2), [-4, 4, -4, 4]);%用法
xlabel('x1');
ylabel('x2');
zlabel('f(x1,x2)');
%非线性规划---料场选址
clc,clear,close all;
a=[1.25,8.75,0.5,5.75,3,7.25];
b=[1.25,0.75,4.75,5,6.5,7.25];
A=[1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0;0 0 0 0 0 0 1 1 1 1 1 1 0 0 0 0];
bb=[20;20];
Aeq=[eye(6),eye(6),zeros(6,4)];
beq=[3;5;4;7;6;11;];
lb=[zeros(12,1);-inf;-inf;-inf;-inf];
x0=[3,5,0,7,0,1,0,0,4,0,6,10,5,1,2,7];
[xans,fvans]=fmincon(@fun2,x0,A,bb,Aeq,beq,lb,[]);