2012年第一届数学建模国际赛小美赛
B题 大规模灭绝尚未到来
原题再现：
亚马逊是地球上现存最大的雨林，比地球上任何地方都有更多的野生动物。它位于南美洲大陆的北侧，共有9个国家：巴西、玻利维亚、厄瓜多尔、秘鲁、哥伦比亚、委内瑞拉、苏里南、法属圭亚那和圭亚那。这类信息需要更加普及，因为砍伐森林正成为森林生存的危险。虽然伐木和森林清理在短期内可能有利可图，但从长远来看，这对亚马逊来说是一种危险。
近年来，森林砍伐率下降到创纪录的低点，巴西50%以上的雨林现在处于某种形式的保护状态。但栖息地丧失的影响需要时间来显现。伦敦帝国理工学院的生态学家罗伯特·埃沃斯说：“砍伐树木并不会直接杀死一只鸟。这些鸟实际上需要很多时间才能死亡。它们都挤进了剩下的栖息地。然后，死亡率会逐渐增加。
要求：
模拟森林砍伐的影响以及在不久的将来是否会发生大规模灭绝。
使用您的模型来估计未来20年内有多少物种将灭绝，请考虑以下四种情况：“一切照旧”，或很少监管治理”，或某种规定根据巴西政府的目标，到2020年将森林砍伐减少80%；到2020年“结束森林砍伐”。
写一份20页的报告（不包括汇总表），介绍你的模型和你对毁林滞后影响的分析。一定要包括政府在确保生态安全方面应该扮演的角色。
整体求解过程概述(摘要)
近年来，在丰富的热带雨林中，特别是在亚马逊河流域，森林砍伐率达到了前所未有的水平，导致生态环境发生变化，许多物种因丧失生活习惯而灭绝。由于亚马逊在巴西的面积最大，本文以巴西的数据为基础进行讨论。到目前为止，已经有许多关于森林砍伐的影响和大灭绝估计的研究。然而，现有的研究大多基于一个方面，为了更深入地探讨这些问题，我们考虑了典型相关分析方法来讨论毁林对经济和生态的影响，并利用偏最小二乘回归方法（PLS）建立了回归函数。
通过求解模型，得到一组典型变量（见方程（1.1）和（1.2）），说明森林砍伐对经济的影响最大，其次是畜牧业。结果表明，当地经济发展的需要，特别是农业经济发展的需要，是热带雨林砍伐的外部驱动因素。基于PLS得到的方程表明，森林砍伐使CO2排放量和濒危物种数量增加，人均可再生淡水资源减少，对生态环境造成严重破坏，对许多野生动物构成巨大威胁。
随后，对于近期是否会发生大灭绝的问题，首先采用非线性函数对毁林曲线进行拟合，拟合函数可以在方程（3.4）中观察到，利用该函数可以预测巴西未来的毁林情况。结合该结果和方程（3.3），预测2014年濒危物种数量为508种，仅比2009年增加1种。考虑到这一点，我们得出结论，在未来可能不会有大规模灭绝。
最后，为了估计未来20年有多少物种将灭绝，我们建立了一个DECS模型来估计未来的灭绝。首先，利用建立的PLS方法计算巴西亚马逊河流域的剩余森林面积。其次，利用微分方程模型计算了时间（t）可能发生的绝灭和“绝灭债”。最后，我们选择计算机模拟的模型来计算在不同的四种情况下灭绝的物种数量。结果表明，BAU、GOV、SR、EOD的灭绝物种分别为13、11、9、7（见表11）。
模型假设：
（1） 亚马逊未来不会发生大火灾。
（2） 未来亚马逊不会发生物种入侵。
（3） 所有数据都是合理的。
（4） 巴西的环境不会受到其他因素的污染。
（5） 这个词只有一个意思。
（6） 巴西的经济不会受到其他国家的严重影响。
问题分析：
由于60～70%的森林被用于养牛，20～30%的森林被用于农业，1～10%的森林被用于伐木[5]。因此，森林砍伐的影响可以通过草原、农业用地和伐木用地的每一种影响来表示。为了衡量森林砍伐对经济的影响，选择了三个指标——GDP、出口和耕地面积。然后得到两组变量。其次，可以考虑典型相关分析来分析毁林对经济的影响，而对于生态影响，则采用PLS方法进行分析。
通过观察森林砍伐随时间的变化，可以直观地看出，近十年来，巴西的生态环境日益恶化。考虑到这一点，我们选择建立一个基于偏最小二乘回归（PLS）的模型来分析毁林与生态环境的相关性。
为了估计未来20年有多少物种将灭绝，我们提出了一种估计未来灭绝的技术。首先，我们考虑使用PSL来计算巴西亚马逊地区的剩余森林面积。其次，利用微分方程和计算机模拟两种方法计算物种灭绝和物种灭绝。
模型的建立与求解整体论文缩略图
全部论文请见下方“ 只会建模 QQ名片” 点击QQ名片即可
部分程序代码：(代码和文档not free)
data ex;
Input x1-x3 y1-y3@@;
cards;
12393.68 4374.24 1458.08 6017 551 5764
24018.06 9825.57 2547.37 5084 582 5796.7
37727.3 16832.18 3482.52 4608 604 5828
52565.94 25865.78 5006.28 5525 731 5900
76734.9 25578.3 8896.8 6638 965 5901.2
88552.32 31253.76 10417.92 8822 1185 5900
91040.67 37572.34 15895.99 10891 1378 5970.8
107750.4 39040 9369.6 13660 1606 5950
104824.02 50721.3 13525.68 16526 1979 6100
112995.2 54732.05 8827.75 15945 1530 6075.3
;
proc cancorr data=ex all;var y1-y3;WITH x1-x3;
run;
clc,clear
format long
pz=[18.226 307.5 42.944 490
36.391 312.5 41.376 494
58.042 319 38.419 494
83.438 321.5 30.68 497
111.21 331.5 29.845 502
130.224 325.5 29.066 502
144.509 338 28.904 504
156.16 368 28.546 505
169.071 357.5 28.159 504
176.555 376 27.932 507];
disp('Mean and standard deviation:');
mu=mean(pz)
sig=std(pz) %Mean and standard deviation
disp('correlation coefficient matrix:');
rr=corrcoef(pz) %correlation coefficient matrix
data=zscore(pz); %data normalization
n=1;m=3; %N is the number of independent variables, m is the number of the
dependent variable
x0=pz(:,1:n);
y0=pz(:,n+1:end); %The original dependent and independent variables and data
e0=data(:,1:n);f0=data(:,n+1:end); %After standardization of dependent and
independent variables and data
num=size(e0,1);%The number of sample points
chg=eye(n); %transformation matrix initialization
for i=1:n
%the score of w，w* and t
matrix=e0'*f0*f0'*e0;
[vec,val]=eig(matrix); %Eigenvalues and eigenvectors
val=diag(val); %Puts forward characteristic value
[val,ind]=sort(val,'descend');
w(:,i)=vec(:,ind(1)); % put forward the corresponding feature vector of
the maximum eigenvalue
w_star(:,i)=chg*w(:,i); %Calculation the value of w *
t(:,i)=e0*w(:,i); %calculate the score of ti
alpha=e0'*t(:,i)/(t(:,i)'*t(:,i)); %calculate the value of alpha_i
chg=chg*(eye(n)-w(:,i)*alpha'); %calculate transformation matrix
e=e0-t(:,i)*alpha'; %Calculating the value of residual matrix
e0=e;
beta=t\f0; %The coefficient of regression equation
cancha=f0-t*beta; %Residual matrix
ss(i)=sum(sum(cancha.^2)); %error sum of squares
%calculate the value of press(i)
for j=1:num
t1=t(:,1:i);f1=f0;
she_t=t1(j,:);she_f=f1(j,:);
t1(j,:)=[];f1(j,:)=[];
beta1=[t1,ones(num-1,1)]\f1;
cancha=she_f-she_t*beta1(1:end-1,:)-beta1(end,:);
press_i(j)=sum(cancha.^2);
end
press(i)=sum(press_i);
Q_h2(1)=1;
if i>1, Q_h2(i)=1-press(i)/ss(i-1); end
if Q_h2(i)<0.0975
fprintf('The number of components is put forward r=%d',i); break
end
end
%--------------------------------------------------------------------------
disp('Y * about t of the regression coefficient:');
beta_z=t\f0
disp('Y * about x of the regression coefficient:');
xishu=w_star*beta_z
mu_x=mu(1:n);mu_y=mu(n+1:end);
sig_x=sig(1:n);sig_y=sig(n+1:end);
ch0=mu_y-(mu_x./sig_x*xishu).*sig_y;
for i=1:m
xish(:,i)=xishu(:,i)./sig_x'*sig_y(i);
end
disp('Display the original data regression equation coefficient:');
sol=[ch0;xish]
