diff --git a/LPej1_3.m b/LPej1_3.m new file mode 100644 index 0000000..2804b6a --- /dev/null +++ b/LPej1_3.m @@ -0,0 +1,40 @@ +clc; clear; close all; format long g; +prob=optimproblem; +x=optimvar('x',9,1,'LowerBound',0,'Type','integer'); +p=[0.25,0.35,0.50]; +r=[1.25,2.00,2.80]; +ot=[5,10,0;7,9,12;6,8,0;4,0,11;7,0,0]; +total=1000:200:3000; +odp=[300/6000,321/10000,250/4000,783/7000,200/4000]; + +X=[]; +Q=[]; + +for i=1:length(total) + % prob.Objective=(ot(1,1)*ones(1,3)*x(1:3)+ot(1,2)*x(7))*odp(1)+(ot(2,1)*ones(1,3)*x(4:6)+ot(2,2)*x(8)+ot(2,3)*x(9))*odp(2)+(ot(3,1)*(x(1)+x(4))+ot(3,2)*(x(7)+x(8)))*odp(3)+(ot(4,1)*(x(2)+x(5))+ot(4,3)*x(9))*odp(4)+(ot(5,1)*(x(3)+x(6)))*odp(5)-(r(1)-p(1))*ones(1,6)*x(1:6)-(r(2)-p(2))*ones(1,2)*x(7:8)-(r(3)-p(3))*x(9); + dp1=(ot(1,1)*ones(1,3)*x(1:3)+ot(1,2)*x(7))*odp(1); + dp2=(ot(2,1)*ones(1,3)*x(4:6)+ot(2,2)*x(8)+ot(2,3)*x(9))*odp(2); + dp3=(ot(3,1)*(x(1)+x(4))+ot(3,2)*(x(7)+x(8)))*odp(3); + dp4=(ot(4,1)*(x(2)+x(5))+ot(4,3)*x(9))*odp(4); + dp5=(ot(5,1)*(x(3)+x(6)))*odp(5); + dp=dp1+dp2+dp3+dp4+dp5; + prob.Objective=dp-(r(1)-p(1))*ones(1,6)*x(1:6)-(r(2)-p(2))*ones(1,2)*x(7:8)-(r(3)-p(3))*x(9); + prob.Constraints.con1=p(1)*ones(1,6)*x(1:6)+p(2)*ones(1,2)*x(7:8)+p(3)*x(9)+dp<=total(i);%总费用 + prob.Constraints.con2=ot(1,1)*ones(1,3)*x(1:3)+ot(1,2)*x(7)<=6000;%A1 + prob.Constraints.con3=ot(2,1)*ones(1,3)*x(4:6)+ot(2,2)*x(8)+ot(2,3)*x(9)<=10000; + prob.Constraints.con4=ot(3,1)*(x(1)+x(4))+ot(3,2)*(x(7)+x(8))<=4000; + prob.Constraints.con5=ot(4,1)*(x(2)+x(5))+ot(4,3)*x(9)<=7000; + prob.Constraints.con6=ot(5,1)*(x(3)+x(6))<=4000; + [sol,fval,flag,out]=solve(prob); + xx=sol.x; + + X=[X,xx]; + Q=[Q,-fval]; + +end +plot(total,Q,'*-'); + + + + + diff --git a/Mengtkl.m b/Mengtkl.m new file mode 100644 index 0000000..f7a9028 --- /dev/null +++ b/Mengtkl.m @@ -0,0 +1,61 @@ + +%27页例题 +clc; clear; close all; +n = 10000; % 使用较小的 n 值以便更容易可视化 +x = unifrnd(0, 12, [1, n]); +y = unifrnd(0, 9, [1, n]); +ans=sum(y < x.^2 & x <= 3)+sum(y < 12 - x & x >= 3); +ans=ans/n; +% 找出满足条件的点 +condition1 = y < x.^2 & x <= 3; +condition2 = y < 12 - x & x >= 3; +condition_met = condition1 | condition2; % 满足任一条件的点 +condition_not_met = ~condition_met; % 不满足任何条件的点 + +% 创建图形窗口 +figure; +hold on;%在同一张图上绘图 + +% 绘制不满足任何条件的点 +scatter(x(condition_not_met), y(condition_not_met), 'k.'); % k----黑色 .----绘制样式 +%scatter绘制散点图 +%x(condition_not_met) 会返回一个新的向量,其中只包含 x 中对应 condition_not_met 为 true 的元素。 + +% 绘制满足第一个条件的点 +scatter(x(condition1), y(condition1), 'r.'); % 红色 + +% 绘制满足第二个条件的点 +scatter(x(condition2), y(condition2), 'b.'); % 蓝色 + +% 添加图例和标签 +legend('不满足任何条件的点', '满足 y < x^2 且 x <= 3 的点', '满足 y < 12 - x 且 x >= 3 的点'); +xlabel('x'); +ylabel('y'); +title('随机生成的点和满足条件的点'); +hold off; + + +%蒙特卡洛法求圆周率qw +clc;clear;close all; + +n=10^5; +x=unifrnd(-1,1,[1,n]); +y=unifrnd(-1,1,[1,n]); +con1=x.^2+y.^2<=1; +con2=~con1; +ans=sum(x.^2+y.^2<=1); +ans=ans/n*4; + +figure ; +hold on; + +scatter(x(con1),y(con1),'r.'); +scatter(x(con2),y(con2),'k.'); + +hold off; + + + + + + diff --git a/NLP.m b/NLP.m new file mode 100644 index 0000000..7f46e9b --- /dev/null +++ b/NLP.m @@ -0,0 +1,135 @@ + +%非线性规划 + +%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=35; +for i=1:5 + con(i+1)=x(i)+x(i+1)>=cnt(i+1); +end +prob.Constraints.con=con; +[sol,fval,flag,out]=solve(prob); +X=sol.x; + + +%linprog解法 +clc; clear; close all; format long g; +f=[1,1,1,1,1,1]; +intcon=[1:6]; +A=zeros(6,6); +A(1,1)=-1; +A(1,6)=-1; +for i=1:5 + A(i+1,i)=-1; + A(i+1,i+1)=-1; +end +lb=zeros(6,1);%注意不是lb=0 +b=[-35;-40;-50;-45;-55;-30]; +[x,fval]=intlinprog(f,intcon,A,b,[],[],lb); + + + +%01规划 +%背包问题 + +f=-[540,200,180,350,60,150,280,450,320,120]; +intcon=1:10; +lb=zeros(10); +ub=ones(10); +A=[6,3,4,5,1,2,3,5,4,2]; +b=30; +[x,fval]=intlinprog(f,intcon,A,b,[],[],lb,ub); + + + + diff --git a/mm1.m b/mm1.m new file mode 100644 index 0000000..fafc5c7 --- /dev/null +++ b/mm1.m @@ -0,0 +1,140 @@ +%投资利益与风险1998A题 + + +%模型一 给定风险承受程度,求最大利益 +f=[-0.05,-0.27,-0.19,-0.185,-0.185]; + +%A矩阵 +%A=[0,0.025,0.015,0.055,0.026]; %错误 +%b=[1,1,1,1,1]; +A=[zeros(4,1),diag([0.025,0.015,0.055,0.026])];%不等式约束条件矩阵 + +%Aeq、beq +Aeq=[1,1.01,1.02,1.045,1.065]; +beq=1; + +%lb +%lb=0; %错误 +lb=zeros(5,1); + + +%可承担风险率a +a=(0:0.001:0.05); + +%保存最优解 +Q=zeros(1,length(a)); +xx=[];%空矩阵存放最优解对应x的值 +for i=1:length(a) + b=a(i)*ones(4,1); + [x,y]=linprog(f,A,b,Aeq,beq,lb); + Q(i)=-y;%注意取负!!! + xx=[xx;x']; +end +plot(a,Q,'*r'); +xlabel("风险率"); +ylabel("最大收益"); + + +%模型二 收益、风险按权重组合 +% f0=[-0.05,-0.27,-0.19,-0.185,-0.185]; +% w=(0:0.1:1); +% Aeq=[1,1.01,1.02,1.045,1.065,0]; +% beq=1; +% lb=0; +% xx=[]; +% Q=zeros(1,length(w)); +% A=[zeros(5,1),diag([0.025,0.025,0.055,0.065,0])]; +% b=ones(5,1); +% for i=1:length(w) +% f=[-w(i)*f0,1-w(i)]; +% b=x(end)*b; +% [x,y]=linprog(f,A,b,Aeq,beq,lb); +% Q(i)=-y; +% xx=[xx,x']; +% end +% plot(w,Q,'*r'); + +%模型二 收益、风险按权重组合 +clc; clear; close all; format long g; + +M = 10000; +prob = optimproblem; +x = optimvar('x', 6, 1, 'LowerBound', 0); +r = [0.05, 0.28, 0.21, 0.23, 0.25]; +p = [0, 0.01, 0.02, 0.045, 0.065]; +q = [0, 0.025, 0.015, 0.055, 0.026]'; +%w = 0:0.1:1; +w = 0.7:0.03:1; +V = []; +Q = []; +X = []; + +prob.Constraints.con1 = (1 + p) * x(1:end-1) == M; +prob.Constraints.con2 = (q(2:end).* x(2:end-1))<= x(end); %下标从1开始 + +for i = 1:length(w) + prob.Objective = w(i) * x(end) - (1 - w(i)) * (r - p) * x(1:end-1); %注意大小写 + [sol, fval, flag, out] = solve(prob); + xx = sol.x; + V = [V, max(q.* xx(1:end-1))]; + Q = [Q, (r - p) * xx(1:end-1)]; + X = [X, xx]; + + plot(V, Q, '*-'); + grid on; + xlabel('风险'); + ylabel('收益'); +end + + + +%模型三:达到一定盈利水平,极小化风险 + +clc; clear; close all; format long g; +M=10000; +k=1500:100:3000; +prob = optimproblem; +x = optimvar('x', 5, 1, 'LowerBound', 0);%下界为0 +r = [0.05, 0.28, 0.21, 0.23, 0.25]; +p = [0, 0.01, 0.02, 0.045, 0.065]; +q = [0, 0.025, 0.015, 0.055, 0.026]'; + +V = []; +Q = []; +X = []; +t = optimvar('t', 'LowerBound', 0); + +%prob.Objective=max(q.*x);%极小化风险 +for i=1:length(k) + prob.Objective = t; % 极小化风险 + prob.Constraints.con1 = (1 + p) * x == M; + prob.Constraints.con2=((r-p)*x>=k(i));%达到一定盈利水平 + prob.Constraints.con3=(q.*x<=t); + [sol,fval,flag,out]=solve(prob); + if flag==1 + xx=sol.x; + X=[X,xx]; + Q=[Q,(r-p)*xx]; + V=[V,fval]; + else + xx=-1*ones(5,1); + X=[X,xx]; + Q=[Q,-1]; + V=[V,-1]; + end + + +end +plot(k,V,'*-'); + + + + + + + + + + + + diff --git a/mm2_graph.m b/mm2_graph.m new file mode 100644 index 0000000..4a0e1a8 --- /dev/null +++ b/mm2_graph.m @@ -0,0 +1,94 @@ +%ej.4.10仓库选址 +clc; clear; close all; format long g; +a=zeros(6); +a(1,[2,5])=[20,15]; +a(2,[1,3,4,5])=[20,20,60,25]; +a(3,[2,4,5])=[20,30,18]; +a(4,[2,3])=[60,30]; +a(5,[1,2,3,6])=[15,25,18,15]; +a(6,5)=15; +s = cellstr(strcat('销售点',int2str([1:6]'))); +G=graph(a,s); +P = plot(G,'layout','force','EdgeColor','k','NodeFontSize',12); +D=[]; +for i=1:6 + temp=0; + for j=1:6 + [path,d]=shortestpath(G,i,j); + if d>temp; + temp=d; + end + end + D=[D,temp]; +end +[minVal,minInd]=min(D); +disp('仓库应建在销售点'); +disp(minInd); + + +%答案版 +clc; clear; close all; + +% 初始化邻接矩阵 +a = zeros(6); +a(1, [2, 5]) = [20, 15]; +a(2, 3:5) = [20, 60, 25]; +a(3, [4, 5]) = [30, 18]; +a(5, 6) = 15; + +% 创建节点标签 +s = cellstr(strcat('V', int2str([1:6]'))); + +% 创建图 +G = graph(a, s, 'upper'); + +% 计算所有节点之间的最短路径距离矩阵 +d = distances(G); + +% 计算从每个节点到其他节点的最大最短路径距离 +D = max(d, [], 2)'; %第三个参数 2 指定沿着矩阵的第二个维度(即列)进行操作。 + +% 找到最大最短路径距离中的最小值及其下标 +[minD, minIndex] = min(D); + +% 绘制图 +plot(G, 'EdgeLabel', G.Edges.Weight, 'Layout', 'force');%设置边的标签 + +% 显示结果 +fprintf('从每个节点到任意其他节点的最大最短路径距离: \n'); +disp(D); +fprintf('最大最短路径距离中的最小值为: %f\n', minD); +fprintf('对应的节点为: V%d\n', minIndex); + +%设备更新问题; +clc;clear;close all; +a=zeros(6); +a(1,[2:6])=[15,20,27,37,54]; +a(2,[3:6])=[15,20,27,37]; +a(3,[4:6])=[16,21,28]; +a(4,[5:6])=[16,21]; +a(5,6)=17; +s=cellstr(strcat(int2str((1:6)'),'年初')); +G=digraph(a,s); +P = plot(G,'layout','force','EdgeColor','k','NodeFontSize',12); +[path,d]=shortestpath(G,1,6);%注意参数个数 +highlight(P,path,"EdgeColor","red",'LineWidth',3.5); + + +%最小生成树 +clc;clear; +a=zeros(9); +a(1,[2:9])=[2 1 3 4 4 2 5 4]; % 顶点1到其他顶点的边的权重 +a(2,[3 9])=[4 1]; % 顶点2到顶点3、顶点9的边的权重 +a(3,4)=1; % 同上。因为写过1到3,和2到3的边的权重,无需重复设 +a(4,5)=1; +a(5,6)=5; +a(6,7)=2; +a(7,8)=3; +a(8,9)=5; +s=cellstr(strcat('节点',int2str([1:9]'))); +G=graph(a,s,'upper'); +P=plot(G,'layout','force'); +T=minspantree(G,'Method','sparse'); +L=sum(T.Edges.Weight); +highlight(P,T,"EdgeColor","red",'LineWidth',2.5);