Compare commits

...

30 Commits

Author SHA1 Message Date
phzrjyvu9 19fa153e5a 因子分析
7 months ago
phzrjyvu9 7c100f95fa 主成分分析
7 months ago
phzrjyvu9 7a93f76da0 层次分析法
7 months ago
p3itgm2rp 5a9d5da45d Merge pull request '和' (#7) from 盘荣博 into main
7 months ago
盘荣博 7f917cf9e6 Merge remote-tracking branch 'remotes/origin/盘荣博' into panrongbo
7 months ago
盘荣博 7c694554d4 修改
7 months ago
p3itgm2rp a1f852b393 Merge pull request '数据可视化' (#6) from 盘荣博 into main
7 months ago
盘荣博 6059412e5a 数据的标准化
7 months ago
盘荣博 e0d39ed553 盘荣博
7 months ago
p3itgm2rp 910206eb1b 添加代码
7 months ago
JCHPJP f12366ca01 acm
7 months ago
p3itgm2rp b19bd9bef5 线性规划的例题
7 months ago
JCHPJP 7b70265a0d 尝试线性模型的acipy模块
7 months ago
JCHPJP ab8bd7dfae 风险投资模型的例子
7 months ago
JCHPJP 22043cdac7 Merge branch '盘荣博' of https://bdgit.educoder.net/p3itgm2rp/mycode
7 months ago
p3itgm2rp 0ec602fbc1 Merge pull request '第一次合并W' (#3) from 王朝群 into main
7 months ago
fighting b789574e3e 第一次尝试
7 months ago
盘荣博 3d0453cf0b 数据可视化
7 months ago
p3itgm2rp aa9c8889bd Merge pull request '三人三个文件夹' (#2) from 盘荣博 into main
7 months ago
盘荣博 19262e6ef2 最小生成树和最短路径
7 months ago
盘荣博 7d1020b93a Initial commit
7 months ago
盘荣博 cb4ffe074c 图(最短路径)(最小生成树)
7 months ago
盘荣博 f5a4c8b16b 投资风险收益、多目标规划(约束法、线性加权)
7 months ago
盘荣博 811bec3b62 整数规划(linprog\optimproblem解法),01规划(背包问题)
7 months ago
盘荣博 d032cbab43 非线性规划料场选址目标函数
7 months ago
盘荣博 97a3f6598d 非线性简单示例目标函数定义
7 months ago
盘荣博 2237326033 非线性规划(P42彩电,求解和灵敏性分析)(简单示例)(选址问题)
7 months ago
盘荣博 2a3a14b6d9 蒙特卡洛(不规则面积和圆周率)
7 months ago
盘荣博 95f484af5c 线性规划例题1_3
7 months ago
p3itgm2rp 371c0c70c3 Merge pull request '第一周' (#1) from 姚安欣 into main
7 months ago

3
.idea/.gitignore vendored

@ -0,0 +1,3 @@
# 默认忽略的文件
/shelf/
/workspace.xml

@ -0,0 +1,12 @@
<component name="InspectionProjectProfileManager">
<profile version="1.0">
<option name="myName" value="Project Default" />
<inspection_tool class="PyPackageRequirementsInspection" enabled="true" level="INFORMATION" enabled_by_default="true">
<option name="ignoredPackages">
<value>
<list size="0" />
</value>
</option>
</inspection_tool>
</profile>
</component>

@ -0,0 +1,6 @@
<component name="InspectionProjectProfileManager">
<settings>
<option name="USE_PROJECT_PROFILE" value="false" />
<version value="1.0" />
</settings>
</component>

@ -0,0 +1,4 @@
<?xml version="1.0" encoding="UTF-8"?>
<project version="4">
<component name="ProjectRootManager" version="2" project-jdk-name="Python 3.11 (mycode)" project-jdk-type="Python SDK" />
</project>

@ -0,0 +1,8 @@
<?xml version="1.0" encoding="UTF-8"?>
<project version="4">
<component name="ProjectModuleManager">
<modules>
<module fileurl="file://$PROJECT_DIR$/.idea/mycode.iml" filepath="$PROJECT_DIR$/.idea/mycode.iml" />
</modules>
</component>
</project>

@ -0,0 +1,8 @@
<?xml version="1.0" encoding="UTF-8"?>
<module type="PYTHON_MODULE" version="4">
<component name="NewModuleRootManager">
<content url="file://$MODULE_DIR$" />
<orderEntry type="inheritedJdk" />
<orderEntry type="sourceFolder" forTests="false" />
</component>
</module>

@ -0,0 +1,6 @@
<?xml version="1.0" encoding="UTF-8"?>
<project version="4">
<component name="VcsDirectoryMappings">
<mapping directory="" vcs="Git" />
</component>
</project>

@ -0,0 +1,2 @@
a=[0]*int(2e4)
print(len(a))

@ -0,0 +1,54 @@
from functools import cmp_to_key
c=[[],[]]
string=[[],[]]
for i in range(2):
a = eval(input())
while a>0:
s=input().split()
s[1],s[2]=eval(s[1]),eval(s[2])
c[i].append(s)
string[i].append(s[0])
a-=1
def cmp(a,b):
if a[1]>b[1]:
return -1
elif a[1]==b[1]:
if a[2]<b[2] :
return -1
else:
return 1
else :
return 1
c[0]=sorted(c[0],key=cmp_to_key(cmp))
# print(c[0])
c[1]=sorted(c[1],key=cmp_to_key(cmp))
# print(c[1])
# m='123'
# n='123'
# print(m==n)
st="lzr010506"
s=set(string[0])&set(string[1])
# for i in c[0]:
# for j in c[1]:
# if i[0] == j[0] :
# s.append(j[0])
# print(s)
ans1=0
ans2=0
i=0
while c[0][i][0] != st:
# print(c[0][i][0])
if c[0][i][0] in s:
ans1+=1
i+=1
ans1= i+1-ans1
i=0
while c[1][i][0]!=st:
# print(c[1][i][0])
if c[1][i][0] in s:
ans2+=1
i+=1
ans2=i+1-ans2
# print(ans1,ans2)
print(min(ans1,ans2))
# print(c[0],'\n',c[1])

@ -0,0 +1,18 @@
import heapq
import ctypes
a=int(input())
b=input().split()
for i in range(len(b)):
b[i]=int(b[i])
heapq.heapify(b)
# print(b)
sum=0
while len(b) != 1:
c,d=heapq.nsmallest(2,b)
# heapq.heappop(b)
# heapq.heappop(b)
b.remove(c)
b.remove(d)
sum+=c+d
heapq.heappush(b,c+d)
print(sum)

@ -0,0 +1,28 @@
import functools
a,b =map(int, input().split(' '))
nums=[]
for i in range(b):
m,n=map(int,input().split(' '))
if m>n:continue;
nums.append([m,n])
def cmp(a,b):#升序
if(a[0]>b[0]):
return 1
else :
return -1#不变
nums=sorted(nums,key=functools.cmp_to_key(cmp))
# for i in nums:
# print(i)
sum = 0
start,end=nums[0][0],nums[0][1]
count=1
while count<len(nums):
if nums[count][0]>=start and nums[count][0]<=end:
if end<=nums[count][1]:
end=nums[count][1]
else :
sum+=end-start+1
start=nums[count][0];end=nums[count][1]
count+=1
sum+=end-start+1
print(a-sum+1)

@ -0,0 +1,7 @@
a,b =map(int, input().split(' '))
nums=[]
for i in range(b):
m,n=map(int,input().split(' '))
nums.extend(range(m,n+1))
nums=set(nums)
print(a+1-len(nums))

@ -0,0 +1,64 @@
%AHP步骤
clc,clear,close all;
A=[1,2,3,5
1/2,1,1/2,2
1/3,2,1,2
1/5,1/2,1/2,1];
[row,col]=size(A);
%判断矩阵一致性检验
n=col;
maxlam=max(eig(A));
RI=[0 0 0.58 0.90 1.12 1.24 1.32 1.41 1.45];
CI=(maxlam-n)/(n-1);
CR=CI/RI(n);
%判断矩阵确定权重
for i=1:col
sumcol=sum(A(:,i));
for j=1:row
A(j,i)=A(j,i)/sumcol;
end
end
weig=zeros(row,1);
for i=1:row
sumrow=sum(A(i,:));
weig(i)=sumrow/n;
end
%各个指标归一化 按列单位化
data= [1686.4 3183 12000 397
903.6 1916.4 3439.6 43
837.6 817.6 4748 1159
824.9 1296.4 12000 442
2110.2 1465.7 6199.5 228];
[rowd,cold]=size(data);
for i=1:cold
sumcold=sum(data(:,i));
for j=1:rowd
data(j,i)=data(j,i)/sumcold;
end
end
%按权重计算分数
score=data*weig;
projectNames={'老番茄','何同学','木鱼水心','凉风','罗翔'};
figure;
bar(score);%条形图
set(gca, 'XTickLabel', projectNames); %每个条形图标签
xlabel('博主');
ylabel('加权总分');
title('得分');
grid on; % 网格线

@ -1,40 +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,'*-');
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,'*-');

@ -1,61 +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;
%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;

@ -1,135 +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<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,[]);
%非线性规划
%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,[]);

@ -0,0 +1,96 @@
clc,clear,format short,close all;
load data_mh.mat;
[n,p]=size(x);
%标准化
X=zscore(x);
%协方差矩阵/相关系数矩阵
R=cov(X);
[V,D]=eig(R);%[特征向量,特征值]
lambda=diag(D);
lambda=lambda(end:-1:1);
total_contri=sum(lambda);
cum_contri=cumsum(lambda);
contri_rate=cum_contri/total_contri;
V1=rot90(V)';
disp(V1);
c1=V1(:,1);
c2=V1(:,2);
XX=X*c1;
YY=X*c2;
figure;
scatter(XX, YY, 'filled');
xlabel('第一主成分');
ylabel('第二主成分');
title('主成分得分图');
grid on;
%散点
figure;
scatter(XX,YY,'k.');
figure;
plot(c1, c2, 'o', 'MarkerSize', 8, 'MarkerFaceColor', 'b');
text(c1, c2, cellstr(num2str((1:p)')), 'VerticalAlignment','bottom', 'HorizontalAlignment','right');
xlabel('第一主成分');
ylabel('第二主成分');
title('主成分载荷图');
grid on;
%主成分得分图Score Plot
X = randn(100, 5);
% 主成分分析
[coeff, score, latent, tsquared, explained, mu] = pca(X);
% 绘制前两个主成分的得分图
figure;
scatter(score(:,1), score(:,2));
xlabel('第一主成分');
ylabel('第二主成分');
title('主成分得分图');
grid on;
% 添加数据点标签
text(score(:,1), score(:,2), num2str((1:size(score,1))'), 'FontSize', 8);

@ -0,0 +1,46 @@
clc,clear,close all;
%相关系数矩阵
r=[ 1.000,0.577,0.509,0.387,0.462
0.577,1.000,1.599,0.389,0.322
0.509,0.599,1.000,0.436,0.426
0.387,0.389,0.436,1.000,0.523
0.462,0.322,0.426,0.523,1.000];
[vec1,val,rate]=pcacov(r);%特征向量、特征值、贡献率
f1=repmat(sign(sum(vec1)),size(vec1,1),1);%调整符号
vec2=vec1.*f1;%是用 .*
f2=repmat(sqrt(val)',size(vec2,1),1);
a=vec2.*f2;%载荷矩阵
a1=a(:,1);
tcha1=diag(r-a1*a1');
a2=a(:,[1,2]);
tcha2=diag(r-a2*a2');
ccha2=r-a2*a2'-diag(tcha2);
con=cumsum(rate);
clc,clear,close all;
load data_mh.mat;
[n,p]=size(x);
%标准化
X=zscore(x);
%相关系数矩阵
r=cov(X);
[vec1,val,rate]=pcacov(r);%特征向量、特征值、贡献率
f1=repmat(sign(sum(vec1)),size(vec1,1),1);%调整符号
vec2=vec1.*f1;%是用 .*
f2=repmat(sqrt(val)',size(vec2,1),1);
a=vec2.*f2;%载荷矩阵
a1=a(:,1);
tcha1=diag(r-a1*a1');
a2=a(:,[1,2]);
tcha2=diag(r-a2*a2');
ccha2=r-a2*a2'-diag(tcha2);
con=cumsum(rate);

@ -1,3 +1,3 @@
function f=fun1(x)
f=x(1)*x(1)+x(2)*x(2)-x(1)*x(2)-2*x(1)-5*x(2);
end
function f=fun1(x)
f=x(1)*x(1)+x(2)*x(2)-x(1)*x(2)-2*x(1)-5*x(2);
end

@ -1,15 +1,15 @@
function f=fun2(x)
a=[1.25,8.75,0.5,5.75,3,7.25];
b=[1.25,0.75,4.75,5,6.5,7.25];
f1 = 0; % 初始化f1
f2 = 0; % 初始化f2
for i=1:6
s=sqrt((a(i)-x(13))^2+(b(i)-x(14))^2);
f1=x(i)*s+f1;
end
for i=7:12
s=sqrt((a(i-6)-x(15))^2+(b(i-6)-x(16))^2);
f2=s*x(i)+f2;
end
f=f1+f2;
function f=fun2(x)
a=[1.25,8.75,0.5,5.75,3,7.25];
b=[1.25,0.75,4.75,5,6.5,7.25];
f1 = 0; % 初始化f1
f2 = 0; % 初始化f2
for i=1:6
s=sqrt((a(i)-x(13))^2+(b(i)-x(14))^2);
f1=x(i)*s+f1;
end
for i=7:12
s=sqrt((a(i-6)-x(15))^2+(b(i-6)-x(16))^2);
f2=s*x(i)+f2;
end
f=f1+f2;
end

@ -1,49 +1,49 @@
%整数规划
%P23 ej2.5
%optimproblem解法
clc; clear; close all; format long g;
prob=optimproblem;
x=optimvar('x',6,'LowerBound',0,'Type','integer');
prob.Objective=sum(x);
cnt=[35,40,50,45,55,30];
con=optimconstr(6);
con(1)=x(1)+x(6)>=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);
%整数规划
%P23 ej2.5
%optimproblem解法
clc; clear; close all; format long g;
prob=optimproblem;
x=optimvar('x',6,'LowerBound',0,'Type','integer');
prob.Objective=sum(x);
cnt=[35,40,50,45,55,30];
con=optimconstr(6);
con(1)=x(1)+x(6)>=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);

@ -1,140 +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,'*-');
%投资利益与风险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,'*-');

@ -1,94 +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);
%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);

@ -0,0 +1,7 @@
**三个文件夹**
数据可视化:暂时没有文件。
图:其中包括最短路径和最小生成树的展示
线性规划: 使用scipy库便且使用scipy的optimize模块下的linprog函数numpy函数统一操作数据以及使用matplotlib库下的pyplot模块

Binary file not shown.

After

Width:  |  Height:  |  Size: 824 KiB

@ -0,0 +1,28 @@
import networkx as nx
import pylab as plt
import numpy as np
p=[25,26,28,31] ;a=[10,14,18,26] ; r=[20,16,13,11]
b= np.zeros((5,5))
for i in range(5):
for j in range(i+1,5):
b[i,j]=p[i]+np.sum(a[0:j-i])-r[j-i-1]
G=nx.DiGraph(b)
print(G)
p=nx.dijkstra_path(G,source=0,target=4,weight='weight')
print("最短路径为:",np.array(p)+1)#下标从零开始
d=nx.dijkstra_path_length(G,source=0,target=4,weight="weight")
print("所求的费用最小值为:",d)
s=dict(zip(range(5),range(1,6)))
plt.rc("font",size=16)
pos=nx.shell_layout((G))#设置布局
print(type(pos),'\npos=',pos)
w=nx.get_edge_attributes(G,"weight")
print(type(w),'\nw=',w)
nx.draw(G,pos,font_weight='bold',labels=s,node_color='r')#绘制点和边
nx.draw_networkx_edge_labels(G,pos,edge_labels=w)#绘制标签
#绘制最短路径
path_edges=list(zip(p,p[1:]))
print(type(path_edges),"\npath_edges=",path_edges)
nx.draw_networkx_edges(G,pos,edgelist=path_edges,edge_color="r",width=1)
plt.savefig("figure10_9.png",dpi=1000);plt.show()

Binary file not shown.

After

Width:  |  Height:  |  Size: 2.7 MiB

@ -0,0 +1,70 @@
import pandas as pd
<<<<<<< HEAD
import numpy as np
from scipy.stats import zscore
from sklearn.decomposition import PCA
=======
from scipy.stats import zscore
>>>>>>> remotes/origin/盘荣博
import matplotlib.pyplot as plt
from matplotlib.pyplot import ylabel
df = pd.read_excel("棉花产量论文作业的数据.xlsx")
# plt.plot(df["年份"],df["单产"])
plt.rcParams['font.sans-serif']="SimHei"
# plt.rcParams['size'] =10
# plt.ylabel('单产')
# plt.xlabel('年份')
# print(df)
d = df.to_numpy()[:,1:]
print(d)
plt.subplot(4,1,1)
plt.scatter(d[:,:1],d[:,1:2],c='r')
ylabel('原始数据'),plt.title("单产和种子费用的关系")
#公式调用标准化,遵守标准正态分布
data = zscore(d)
print(data)
plt.subplot(4,1,2)
plt.scatter(data[:,:1],data[:,1:2],c='b',)
ylabel('zscore')
print(d.max(axis=0))
print(d.std(axis=0))
print(d.mean(axis=0))
#手写标准正态分布
data1=(d-d.mean(axis=0))/d.std(axis=0)
print(data1)
plt.subplot(4,1,3)
plt.scatter(data1[:,:1],data1[:,1:2],c='y')
ylabel('手写标准正态分布')
data2=(d-d.min(axis=0))/(d.max(axis=0)-d.min(axis=0))
plt.subplot(4,1,4)
plt.scatter(data2[:,:1],data2[:,1:2],c='g')
plt.xlabel('压缩到0~1')
print(data==data1)
<<<<<<< HEAD
# plt.savefig("shuju.jpg",dpi=2000)
# plt.show()
md= PCA().fit(data)
cf = np.cov(data.T)#求协方差矩阵
print(cf)
c, d= np.linalg.eig(cf)
print("特征值:\n",c)
print(md.explained_variance_)
e=c/c.sum()
# for _ in range(len(e)):
# if(_!=0):
# e[_]+=e[_-1]
print('贡献率:')
print(e)
print(md.explained_variance_ratio_)
print('特征向量:')
print(d.T)
print(md.components_)
print(md.components_-d.T<=0.1)
=======
plt.savefig("shuju.jpg",dpi=2000)
plt.show()
>>>>>>> remotes/origin/盘荣博

@ -0,0 +1,19 @@
from matplotlib import pyplot as plt
from numpy import ones , diag , c_ , zeros
from scipy.optimize import linprog
import time
start = time.time()
c=list( [-0.05,-0.27,-0.19,-0.185,-0.185])
A=c_[zeros(4),diag([0.025,0.015,0.055,0.026])]
Aeq = [[1,1.01,1.02,1.045,1.065]];beq=[[1]]
a=0;aa=[];ss = []
while a<=0.05:
b=list(ones(4)*a)
res= linprog(c,A,b,Aeq,beq)
aa.append(a);ss.append(-res.fun)
a=a+0.001
end = time.time()
print("花费时间:",end-start)
plt.plot(aa,ss,"r*")
plt.xlabel("$a$");plt.ylabel("$Q$",rotation=90)
plt.savefig("figure5_1_1.png",dpi=500) ;plt.show()

@ -0,0 +1,26 @@
import numpy as np
from numpy import ones , zeros , c_,diag
from scipy.optimize import linprog
import matplotlib.pyplot as plt
c = np.append(zeros(5).tolist(),[1]).tolist()
print(c)
A=np.append(zeros(4).reshape(4,1),diag([0.025,0.015,0.055,0.026]),axis=1)
A=np.append(A,ones(4).reshape(4,1)*-1,axis=1).tolist()
Aeq =[[1,1.01,1.02,1.045,1.065,0]] ;beq=[1]
A.append([-0.05,-0.27,-0.19,0.185,-0.185,0])
print(A)
k=0.05;step = 0.005
b=([0]*4);b.append(-k)
print(b)
kk=[];ss=[]
while k<0.28:
res= linprog(c,A,b,Aeq,beq)
kk.append(k)
ss.append(res.fun)
print(res.fun)
k+=step
b[4]=-k
plt.plot(kk,ss,'r*')
plt.xlabel("$k$");plt.ylabel('$R$')
plt.savefig("figures5_1_2.png",dpi=500);plt.show()

Binary file not shown.

After

Width:  |  Height:  |  Size: 93 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 115 KiB

@ -0,0 +1,56 @@
from scipy.optimize import linprog
import time
import numpy as np
'''
目标函数
min z= -x1 + 4 x2
约束条件
-3x1 + x2 < = 6
x1 + 2x2 <= 4
x2 >= -3
'''
begin = time.time()
c=[-1,4]
A=[[-3,1],[1,2]]
b=[[6],[4]]
bounds=((None , None),(-3,None))
res= linprog(c,A,b,None,None,bounds)
end= time.time()
print("结果为:",res.fun)
print("x的值位",res.x," ",type(res.x))
print(end-begin,"\n",'----------------------------------------')
'''
线性规划问题
max z = x1 - 2x2 - 3x3
约束条件
-2x1 + x2 + x3 <= 9
-3x1 + x2 + 2x3 >= 4
4x1 -2x2 -x3 = -6
x1 >= -10 ,x2 >=0
'''
'''
scipy标准形式
min w= -x1 + 2x2 + 3x3
约束
-2x1 + x2 + x3 <= 9
3x1 - x2 - 2x3 <= -4
'''
c= (-np.array([1,-2,-3])).tolist()
A=[[-2,1,1],[-3,1,2]]
A[1]= (-np.array(A)[1]).tolist()
b= [[9],[4]]
b[1]=(-np.array(b[1])).tolist()
Aeq=[[4,-2,-1]]
beq =[[-6]]
LB =[-10,0,None]
UB = [None]*len(c)
bounds= tuple(zip(LB,UB))
res= linprog(c,A,b,Aeq,beq,bounds)
end= time.time()
print("结果为:",res.fun)
print("x的值位",res.x," ",type(res.x))
print(end-begin,"\n",'----------------------------------------')
print(res)

@ -0,0 +1,23 @@
import time
import numpy as np
from scipy.optimize import linprog
start = time.time()
c = [-110,-120,-130,-110,-115,150]
c = (-np.array(c)).tolist()
A=[[1,1,0,0,0,0],
[0,0,1,1,1,0],
[8.8,6.1,2.0,4.2,5.0,-6],
[8.8,6.1,2.0,4.2,5.0,-3]
]
A[3]=(-np.array(A[3])).tolist()
b=[[200],[250],[0],[0]]
Aeq =[[1,1,1,1,1,-1]]
beq = [[0]]
LB= [0]*len(c)
UB= [None]*len(c)
bounds= tuple(zip(LB,UB))
res = linprog(c,A,b,Aeq,beq,bounds)
end = time.time()
print("最优解:\n",res.x)
print("目标函数最小值:\n",-res.fun)

@ -0,0 +1,16 @@
import numpy as np
import pylab as pl
from scipy import interpolate
import matplotlib.pyplot as plt
x = np.linspace(0,2*np.pi+np.pi/4,10)
x1 = np.linspace(0,2*np.pi+np.pi/4,100)#num个0~2*Pi+Pi/4的范围的点
y = np.sin(x)
y1= np.sin(x1)
plt.xlabel(f"安培/A")
plt.ylabel(f'伏特/V')
linear_ = interpolate.interp1d(x,y)
print(linear_(x1))
plt.rcParams['font.sans-serif'] = ['SimSun']#设置字体
plt.plot(x,y,'o',label=f"原始数据")
plt.plot(x1,linear_(x1),"*",label="线性插入")
pl.show()

@ -0,0 +1,27 @@
from scipy import optimize
import numpy as np
c= np.array([-1,4])
A=np.array([[-3,1],[1,2]])
b=np.array([6,4])#小于关系
Aeq=np.array([[1,1,1]])#相等
beq = np.array([7])#
bounds = ((None,None),(-3,None))
res = optimize.linprog(c,A,b,None,None,bounds=bounds)
print("目标函数最小值:",res.fun)#目标函数最优解
print("最优解",res.x)#求得的最优解
c = [-1,2,3]
A=[[-2,1,1],[3,-1,-2]]
b=[[9],[-4]]
Aeq =[[4,-2,-1]]
beq=[-6]
LB=[-10,0,None]
UB = [None]*len(c)
print(UB)
bound = tuple(zip(LB,UB))
print(zip(LB,UB),'\n',bound)
res = optimize.linprog(c,A,b,Aeq,beq,bounds=bound)
print("函数的最小值为",res.fun)
print("最优解为:",res.x)
Loading…
Cancel
Save