2023年国赛高教杯数学建模
D题 圈养湖羊的空间利用率
原题再现
规模化的圈养养殖场通常根据牲畜的性别和生长阶段分群饲养，适应不同种类、不同阶段的牲畜对空间的不同要求，以保障牲畜安全和健康；与此同时，也要尽量减少空间闲置所造成的资源浪费。在实际运营中，还需要考虑市场上饲料价格和产品销售价格的波动以及气候、疾病、种畜淘汰、更新等诸多复杂且关联的因素，但空间利用率是相对独立并影响养殖场经营效益的重要问题。
湖羊是国家级绵羊保护品种，具有早期生长快、性成熟早、四季发情并且可以圈养等优良特性。湖羊养殖场通常建有若干标准羊栏，每一标准羊栏所能容纳的羊只数量由羊的性别、大小、生长阶段决定。
湖羊养殖的生产过程主要包括繁殖和育肥两大环节。人工授精技术要求高，因此湖羊繁殖大多采用种公羊和基础母羊自然交配的方式。怀孕母羊分娩后给羔羊哺乳，羔羊断奶后独立喂饲，育肥长成后出栏。自然交配时将若干基础母羊与一只种公羊关在一个羊栏中，自然交配期约为 3 周，然后将种公羊移出。受孕母羊的孕期约为 5 个月，每胎通常产羔 2 只。母羊分娩后哺乳期通常控制在 6 周左右，断奶后将羔羊移至育肥羊栏喂饲。一般情况下，羔羊断奶后经过7 个月左右育肥就可以出栏。母羊停止哺乳后，经过约 3 周的空怀休整期，一般会很快发情，可以再次配种。按上述周期，正常情况下，每只基础母羊每 2 年可生产 3 胎。在不考虑种公羊配种能力差异的情况下，种公羊与基础母羊一般按不低于 1:50 的比例配置。种公羊和母羊在非交配期原则上不关在同一栏中。
某湖羊养殖场设置标准羊栏，规格是：空怀休整期每栏基础母羊不超过 14 只；非交配期的种公羊每栏不超过 4 只；自然交配期每栏 1 只种公羊及不超过 14 只基础母羊；怀孕期每栏不超过 8 只待产母羊；分娩后的哺乳期，每栏不超过 6 只母羊及它们的羔羊；育肥期每栏不超过 14只羔羊。原则上不同阶段的羊只不能同栏。
养殖场的经营管理者为保障效益，需要通过制定生产计划来优化养殖场的空间利用率。这里的生产计划，主要是决定什么时间开始对多少可配种的基础母羊进行配种，控制羊只的繁育期，进而调节对羊栏的需求量，以确保有足够多的羊栏，同时尽量减少羊栏闲置。当羊栏不够时，可以租用其他场地。
请建立数学模型讨论并解决以下问题：
问题 1 不考虑不确定因素和种羊的淘汰更新，假定自然交配期 20 天，母羊都能受孕，孕期 149 天，每胎产羔 2 只，哺乳期 40 天，羔羊育肥期 210 天，母羊空怀休整期 20 天。该湖羊养殖场现有 112 个标准羊栏，在实现连续生产的条件下，试确定养殖场种公羊与基础母羊的合理数量，并估算年化出栏羊只数量的范围。若该养殖场希望每年出栏不少于 1500 只羊，试估算现有标准羊栏数量的缺口。
问题 2 在问题 1 的基础上，对 112 个标准羊栏给出具体的生产计划（包括种公羊与基础母羊的配种时机和数量、羊栏的使用方案、年化出栏羊只数量等），使得年化出栏羊只数量最大。
问题 3 问题 1 和问题 2 中用到的数据都没有考虑不确定性，一旦决定了什么时间开始对多少可配种的基础母羊进行配种，后续对羊栏的安排和需求也就随之确定。例如，用 3 个羊栏给 42 只母羊进行配种，孕期需要 6 个羊栏，哺乳期需要 7 个羊栏给怀孕母羊分娩和哺乳，哺乳期结束就需要给 84 只断奶羔羊和 42 只母羊共安排 9 个羊栏进行育肥和休整。但实际情况并非如此，配种成功率、分娩羔羊的数目和死亡率等都有不确定性，哺乳时间也可以调控，这些都会影响空间需求。
现根据经验作以下考虑：
(1) 母羊通过自然交配受孕率为 85%，交配期结束后 30 天可识别出是否成功受孕；
(2) 在自然交配的 20 天中受孕母羊的受孕时间并不确知，而孕期会在 147-150 天内波动，这些因素将影响到预产期范围；
(3) 怀孕母羊分娩时一般每胎产羔 2 只，少部分每胎产羔 1 只或 3 只及以上，目前尚没有实用手段控制或提前得知产羔数。羔羊出生时，有夭折的可能，多羔死亡率高于正常。通常可以按平均每胎产羔 2.2 只、羔羊平均死亡率 3%估算。
(4) 母羊哺乳期过短不利于羔羊后期的生长，通常是羔羊体重达到一定标准后断奶；而哺乳期过长，母羊的身体消耗就越大，早点断奶，有利于早恢复、早发情配种。一种经验做法是将哺乳期控制在 35-45 天内，以 40 天为基准，哺乳期每减少 1 天，羔羊的育肥期增加 2 天；哺乳期每增加 1 天，羔羊的育肥期减少 2 天。除此之外，母羊的空怀休整期可在不少于 18 天的前提下灵活调控。
此外，如有必要，允许分娩日期相差不超过 7 天的哺乳期母羊及所产羔羊同栏，允许断奶日期相差不超过 7 天的育肥期羔羊同栏，允许断奶日期相差不超过 7 天的休整期母羊同栏。为简化问题，不考虑母羊流产、死亡以及羔羊在哺乳期或育肥期夭折和个体发育快慢等情况。
在以上不确定性的考虑下，生产计划的制定与问题 1 和问题 2 将有较大的不同：一旦作出了“什么时间开始对多少可配种的基础母羊进行配种”的决定，后续羊栏的需求和安排不再是随之确定的，而是每一步都会出现若干种可能的情况需要作相应的并遵从基本规则的安排处理，但无法改变或调整上一步。因此，某种意义上，本问题要讨论研究的生产计划将是一个应对多种可能情况的“预案集”。
请综合考虑可行性和年化出栏羊只数量，制定具体的生产计划，使得整体方案的期望损失最小。其中整体方案的损失由羊栏使用情况决定，当羊栏空置时，每栏每天的损失为 1；当羊栏数量不够时，所缺的羊栏每栏每天的损失（即租用费）为 3。
整体求解过程概述(摘要)
湖羊是优秀的养殖品种，湖羊养殖场通常有若干标准羊栏，根据湖羊的性别和生长阶段分群饲养，不同阶段、性别和大小的羊只对空间要求不同，所以每一羊栏所容纳的羊只数量由上述因素决定，在实际运营中，空间利用率是相对独立并影响养殖场经营效益的重要问题，本文主要研究了湖羊养殖过程中空间利用率的问题，给出了具体生产计划，较好地解决了提出的三个问题。
针对问题 1，母羊的工作周期包括交配期、孕期、哺乳期和休整期，一整个周期持续229天，仅考虑单个批次一定数量的基础母羊以固定周期的方式重复工作，可求出其(包括羔羊)所占用的羊栏天数，在理想状态下多批次交替工作可以把羊栏使用数的波动抹平，而出栏羊数正比于基础母羊的数量n，可求出n对应的每一栏能转化为多少年化出栏数.遍历n 便可以得到年化出栏羊只数量范围的估计为1163 至1312，要想年化出栏羊只数量达到1500，缺口约为16至32。
针对问题2，我们考虑了一般情形，即决策变量为批次之间的间隔g_i和每批次进入交配期的基础母羊数量x_i，以母羊总数量为目标函数，112 个羊栏数量为约束条件建立规划模型，然而该模型是非线性的整数规划，且非线性约束条件非常多，基本不可解。我们对模型进行了化简，固定间隔和数量，即g_i=g，x_i=x，这样决策变量只有两个整数，决策空间有限，用遍历的方法得到最优解g=22,x=40 在一个工作周期内重复 10次，该方案的年化出栏羊只数量1200(2年3胎)。然而该方案最多只使用了110个羊栏，有2个羊栏的元余,且空间利用率只有 95.56%我们经过一定范围的遍历得到了更优解g=22，x_1=48，x_2=…=x_10=40，年化出栏羊只数量达到了1224只,且空间利用率97.38%.为最优解。
针对问题3，我们研究了随机因素对母羊和羔羊各个时期的影响，怀孕时间和孕期的分布使得我们需要在孕期、哺乳期、休整期对母羊分批次管理，包括育肥期的羔羊，不同分支的决策各有不同，我们确定225天的大周期，且每批次之间间隔 25 天，随着时间发展将出现分支，第 51 天，部分没有成功怀孕的母羊退出本次工作，随后按孕期结束时间不同分了3个分支，每个分支尽量保证哺乳期的时长，但不能过于压缩休整期与之对应的羔羊也分为了3个分支，这样只需要确定每批次进入交配期的基础母羊的数量 x，具体方案就确定下来了，我们用计算机模拟充分多的周期，用蒙特卡洛方法计算损失数期望，并以其最小为优化目标，建立优化问题模型，由于决策变量是1维的，我们用遍历法求解我们观察到随着 x的增加损失先递减后递增，在x=40时得到最小的日均损失数3.79，并给出了局部和全局的羊栏数使用情况可视化结果。
模型假设：
1、简单起见，在方案中各批次之间的交配期不重合。
2、按2年3胎的方式计算年化出栏羊只数量。
3、假设羔羊进入育肥期就不会死亡。
4、允许同批次的母羊之间的移动，即从1个该批次的羊栏移至另 1个同批次的羊栏，这样如果同批次的不同羊栏都有母羊移除至下一个阶段，那么剩余的同批次羊可以合并。
5、母羊的怀孕时间和孕期近似服从均匀分布。
问题分析：
问题 1的分析
针对问题 1，注意到母羊的工作周期为 229 天，我们考虑让n 只基础母羊进入交配期，且过了休整期后立刻进入下一个工作周期.在这一个周期内这n只母羊与其生产的羔羊在每个时间段使用多少羊栏是唯一确定的.如果只有1个批次反复来回，那么使用的羊栏数必然波动较大，增加批次可以使得羊栏的使用更平均，我们按理想状态估算即多批次交替后，综合羊栏使用数量关于时间成近似的常数函数关系，那么该批次平均占用的羊栏数量就可以求出来，而年化出栏数与母羊数量正相关，因此可以求出在不同决策下每一个羊栏对应的出栏羊只数量，让 n 从一个固定的范围遍历便可以得到112 个羊栏的条件下出栏羊只数量的范围，也可以估算要达到 1500 只年化出栏羊数额外需要的羊栏数。
问题2的分析
针对问题 2，沿用问题1中固定数量的基础母羊和批次之间固定间隔的生产模式而当每批次基础母羊和间隔给定后，我们需要计算在这样的决策下进入稳定期后，每一天所使用的羊栏数，我们建立优化模型，以母羊数量总和为优化目标，以最大羊栏数为约束条件，而经过简化，决策变量就是二维的整数，我们考虑采用遍历的方法进行求解，并进一步的对我们的结果进行优化。
问题3的分析
针对问题3，仍然采用问题1和2的固定数量固定间隔的生产模式，但是因为有随机因素的参与，问题变得更加复杂。虽然交配期没变，但是孕期结束的时间并不固定,同时还有部分未能成功受孕的母羊因为要进入后续批次而退出，根据我们的简化假设孕期结束的时间相互相差 21 天，而根据孕期哺乳期的条件，需要把同一批进入交配期的羊分成3个分支，这3个分支的状态和决策各不相同，包括后续产下的羔羊也对应的分为 3 个分支。
经过我们分析，应尽可能的延长哺乳期的时间，但是又不能压缩休整期的时间.多方考虑下我们需要确定好工作周期的长度和固定间隔，当基础母羊数量给定后，对应的羊栏使用数量以及相应的损失也可以求出，类似问题 1，让母羊数量遍历找到使得损失最小的方案以最大程度的利用空间。
模型的建立与求解整体论文缩略图
全部论文请见下方“ 只会建模 QQ名片” 点击QQ名片即可
程序代码：(代码和文档not free)
function yi = rdnump(xi,p,y)
yi(1:length(xi))= 0;
pp = yi;
for i= 1:length(p)
pp(i)= sum(p(1:i));
end
for i = 1:y
rr = rand();
index = find(pp>rr,1);
yi(index) = yi(index) + 1;
end
end
clear;
pp=4;
g(1:10)= 22;
p=[20，149，40，2e];
k=10; % floor(sum(p)/p(1)) - 1;
p5 = 21;
x(1:k)= 40;x(1)= 48;
m= ceil(sum(x) / 5e);
m2= ceil(m/4);
nn =[14，8，6，14];
t(1:k,1:4*pp +1)= 0;
T(1:k,1:2*pp)=0;
t(1,1)=1;
forj= 2:5
t(1,j)= t(1,]-1)+ p(j - 1);
end
T(1,1)= t(1,4);T(1,2)= t(1,4) + p5;
for i= 2:pp
t(1,4*(i- 1)+2:4*1+1) = t(1,4*(i - 2)+2:4*(i - 1)+1)+ Sum(p);
T(1,2*(i-1)+1:2*)= T(,2*(i-2)+1:2*(- 1))+ Sum(p);
end
tmax = 1000;
for i=2:k
t(i,:)=t(i-1,:)+g(i-1);
T(i,:)= T(i- 1,:)+ g(i-1);
end
deltaij(1:k,1:4,1:tmax)= 0;
for 1 = 1:tmax
for i = 1:k
kk = find(t(i,:) <= l,1,"last");
j=1+ mod((kk -1),4);
deltaij(i,j,1) =1;
end
end
for 1= 1:sum(p)
for i = 1:k
if max(deltaij(i,:,1))== 0
deltaij(i,41)=1;
end
end
end
deltai bar(1:k,1:tmax)=0;
for l = 1:tmax
for i=1:k
kk = find(T(i,:)<= l,1,"last");
if mod(kk,2) == 1
deltai_bar(i,1) =1;
end
end
end
Delta(1:tmax)=0;
for 1 = 1:tmax
Delta(1) = max(deltaij(:,1,1));
end
N(1:tmax) = 0 ;
for l = 1:tmax
N1=0;N2=0;
for i=1:k
[n1,m1]=numjp(x(i),m);
n=[n1 + ml,ceil(x(i) / nn(2)),ceil(x(i) / nn(3)),ceil(x(i) / nn(4))];
for j =1:4
N1 = N1 + ceil(deltaij(i,j,1)* n(j));
end
N2 = N2 + ceil(2* deltai_bar(i,1) * x(i) / 14);
end
N3=(1 - Delta(1)) * m2;
N(1)=N1+N2+N3;
end
plot((1:length(N)),N)
xlabel('天数')
ylabel('羊栏数')
N max = max(N);
index1 = find(N == max(N),1);
pct = sum(N(indexl:index1 + sum(p) - 1)) / (sum(p)* N max);
NN1(1:k,1:sum(p))= 0;
NN2(1:k,1:sum(p))= 0;
NN3(1:sum(p))= 0;
for l = index1:index1 + sum(p) - 1
for i = 1:k
[n1,m1]= numjp(x(i),m);
n=[n1,ceil(x(i) / nn(2)),ceil(x(i) / nn(3))，ceil(x(i) / nn(4))];
for j = 1:4
NN1(i,l - index1 +1)= NN1(i,l - index1 + 1) + deltaij(i,j,l) * n(j);
end
NN2(i,l - index1 + 1)= NN2(,l - ndex1 + 1) + deltai bar(i,l) * ceil(2 * x(i)
NN3(I index1 +1)= NN3( - index1 + 1) + deltaij(i,1,l)* m1;
end
NN3(1 - index1 +1)= NN3(1 - index1+1)+(1- Delta(1))* m2;
/14);
end
NN_re(1:2*k+1,1:sum(p))= 0 ;
NN_re(1:2:2*k,:)= NN1;
NN re(2:2:2*k,:)= NN2;
NN re(2* k +1,:)= NN3;
clear;
c_f=[];
Gf=[];
for x= 30:60
P= 1000;
imax=9 * p - 1;
T=25*(9*P-1)+ 42;
N(1:1:T)=0;
C= N;
Deltax=[];
G_bar =[];
i=1;
[NN,Delta xx,GG bar] = f3(x);
Delta x =[Delta x,Delta xx];
G bar =[G bar,GG bar];
N(1+25*(i-1):14 25*(i - 1)+42)= N(1+25*(i- 1):1+25*(i-1)420) + NN;
i=2;
[NN,Delta xx,GG bar] = f3(x)
Delta x =[Delta x,Delta xx];
G bar = [G bar,GG bar];
N(1+25*(i- 1):1+25* (i - 1)+42)= N(1+25*(i- 1):1+25* (i- 1)428)+NN;
for i = 3:imax
xx=x + Delta x(i-2);
[NN,Delta xx,GG bar] = f3(xx);
Delta x =[Delta x,Delta xx];
G bar =[G bar,G bar];
N(1+25*(i- 1):1+25*(i-1)+42) = N(1+ 25*(i-1):1+25*(i-1)+ 42)+NN;
end
index1 = 901:1808;
index2=451:225*P;
for t = 1:T
if N(t) >=112
C(t)=3*(N(t)-112);
else
c(t)=112 - N(t);
end
end
c_f =[C_f;mean(C(index2))];
[G_f;sum(G_bar(18:9*P-1)) /(6*P-12)];
end
plot(index1,N(index1))
xlabel('天数')
ylabel('羊栏数')
