一、填空題
- 在種群增長預測問題中,若資源環境等因素是有限的,則應使用的微分方程模型為 Logistic模型
- 某種群分為 4 個年齡組, 各組的繁殖率分別為 0, 0.8, 1.8, 0.2, 存活率分別為 0.5, 0.7, 0.9, 0. 現各組的數量均為 100, 則該種群的的穩定分布向量為
解:穩定分布向量為矩陣最大特征根所對應的特征向量歸一化以后的向量【數學模型P203-206】,下面用R求解
> # --第二題--
> # 構造矩陣L
> L <- matrix(0, ncol = 4, nrow = 4)
> L[1, ] <- c(0, 0.8, 1.8, 0.2)
> L[2, 1] = 0.5
> L[3, 2] = 0.7
> L[4, 3] = 0.9
> eigen(L) -> L.eig
> Re(L.eig$vec[, 1]/sum(L.eig$vec[, 1])) -> w
> w
[1] 0.4769818 0.2307251 0.1562485 0.1360446
于是穩定分布向量為 \((0.4769818, 0.2307251, 0.1562485, 0.1360446)\).
- 工藝品店前段時間的銷售額為
19 10 9 11 9 15 12 17 12 10 11 9 12 8 10 11 15 11 10 14 11 15 10 6 16 15 11 16 13 12 11
已知每個工藝品的成本是 300 元,售價是 800 元,但在關店時仍未售出則只能丟棄.若希望收益最大,那么每日制作工藝品數量應為_______.
解:【數學模型P282】
將工藝品的需求視為連續型隨機變數,用概率密度函式描述概率分布
> shapiro.test(xse)
Shapiro-Wilk normality test
data: xse
W = 0.95367, p-value = https://www.cnblogs.com/hznudmh/p/0.197
shapiro檢驗結果可以認為需求量服從正態分布. 設其平均值為 \(\mu\), 標準差為 \(\sigma\),則需求量為 \(N(\mu,\sigma^2)\)
設每日制作的工業品數量為 \(q\),\(s_1\) 為每日售出一單位工藝品的獲利,\(s_2\) 為因未售出而丟棄的一單位工藝品的損失,\(r\) 為某天的需求量. 記需求量的概率密度函式為 \(p(r)\). 于是日均利潤為
\[E(q) = \int_{-\infty}^q[s_1r-s_2(q-r)]p(r)\mathrm{d}r + \int_{q}^{\infty}s_1qp(r)\mathrm{d}r \]則要求 \(E(q)\) 最大,令 \(\frac{\mathrm{d}E}{\mathrm{d}q}=0\),記 \(r\) 的概率分布函式為 \(F(q) = \int_{-\infty}^qp(r)\mathrm{d}r\),可得 \(q\) 的最優解滿足
\[F(q) = \frac{s_1}{s_1+s_2} \]根據題意 \(s_1 = 800-300=500,\ s_2 = 300\). 下面用R求q的值
> s1 <- 500
> s2 <- 300
> mean(xse) -> mu
> sd(xse) -> sigma
> f <- s1/(s1+s2)
> qnorm(f, mu, sigma)
[1] 12.89302
故每日制作的最佳工藝品數量為13
二、建模與計算題
- 某陣地(原點)向位于(L,0)處的敵艦發射導彈,敵艦隨即以速度(0,v)逃竄,若在任何時刻(位置),導彈的飛行方向都對準敵艦,且速度大于 v,請列出其飛行曲線 y(x)所滿足的方
程.
解:記 \(x = x(t),\ y = y(t)\),設導彈合速度為 \(u\),將 \(u\) 分解為水平方向的 \(u_x\) 和豎直方向的 \(u_y\). 示意圖如下(有點丑別介意)
容易列出
- 每隔 1 小時測得的某細菌的生物量資料為:9.6, 18.3, 29.0, 47.2, 71.1, 119.1, 174.6, 257.3, 350.7, 441.0, 513.3, 559.7, 594.8, 629.4, 640.8, 651.1, 655.9, 659.6, 661.8. 試建立 Logistic 模型,并給出下 1 個小時該細菌生物量的預測值
解. 設生物量為 \(x\),生物量增長率為 \(r(x)\),不妨設 \(r(x) = a + bx\). 記 \(x_m\) 為所能容納的生物量最大值,且當 \(x = x_m\) 時細菌量不再增長;\(r\) 為 \(x=0\) 時的增長率.
于是可以得到 \(a = r, \ b = -\frac{r}{x_m}\). 可以得到 Logistic 模型
那么只需估計出引數 \(x_m\) 和 \(r\) 即可.
- 方法一 線性最小二乘法估計
方程式改寫為
R代碼:
> # --第二題--
> x <- c(9.6, 18.3, 29.0, 47.2, 71.1, 119.1,
+ 174.6, 257.3, 350.7, 441.0, 513.3, 559.7,
+ 594.8, 629.4, 640.8, 651.1, 655.9, 659.6, 661.8)
> # 線性最小二乘估計
> # 使用數值微分近似估計增長率
> dx <- c()
> dx[1] <- (-3*x[1]+4*x[2]-x[3])/2
> for (i in 2:(length(x)-1)){
+ dx[i] <- (x[i+1]-x[i-1])/2
+ }
> dx[length(x)] <- (x[length(x)-2]-4*x[length(x)-1]+3*x[length(x)])/2
> xx <- dx/x
> # 最小二乘估計
> lm(xx~x) -> sol2
> summary(sol2)
Call:
lm(formula = xx ~ x)
Residuals:
Min 1Q Median 3Q Max
-0.088687 -0.019463 -0.006286 0.006698 0.234409
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.761e-01 2.579e-02 22.34 4.88e-14 ***
x -8.780e-04 5.692e-05 -15.43 1.98e-11 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.06389 on 17 degrees of freedom
Multiple R-squared: 0.9333, Adjusted R-squared: 0.9294
F-statistic: 238 on 1 and 17 DF, p-value: 1.983e-11
> 5.761e-01/8.780e-04
[1] 656.1503
于是 \(r = 0.5761, \ x_m = 656.1503\). 下面用 R 預測
> xm <- 5.761e-01/8.780e-04
> x0 <- 9.6
> r <- 5.761e-01
> f = function(t) xm/(1+(xm/x0-1)*exp(-r*t))
> f((length(x)+1))
[1] 655.7127
故下一個小時的預測值為655.7127.
- 方法二 非線性最小二乘估計
> library(pracma)
> Ffun<-function(y){
+ f<-numeric(length(x))
+ for(i in 1:(length(x))) {
+ f[i]<-x[i]-y[2]/(1+(y[2]/y[1]-1)*exp(-y[3]*(i-1)))
+ }
+ f
+ }
> sol22<-lsqnonlin(Ffun,x0=c(10,666,0.57))
> sol22
$x
[1] 9.1355173 663.0219910 0.5469947
$ssq
[1] 194.3254
$ng
[1] 0.0007158751
$nh
[1] 3.155265e-07
$mu
[1] 0.004704803
$neval
[1] 15
$errno
[1] 1
$errmess
[1] "Stopped by small x-step."
于是 \(r = 0.5469947,\ x_m = 663.0219910\).
預測:
> x0 <- sol22$x[1]
> xm <- sol22$x[2]
> r <- sol22$x[3]
> f(length(x)+1)
[1] 662.1814
故下一個小時的預測值為662.1814.
- 某旅行團計劃進行一次假期旅行, 有 4 個備選方案. 已知選擇目的地時的參考指標有費
用(萬元),景色(分),食宿條件(分),勞累程度(分)四項, 4 處目的地的得分如下表所示. 請使用理想數法選擇最適合的目的地.(矩陣標準化時指定使用模一化方法)
解:先構造決策矩陣,并模一化
由于“費用”和“勞累程度”為費用型屬性,因此要取倒數化為效益型屬性
R代碼:
> # --第三題--
> D <- matrix(c(1/5000, 1/3000, 1/4000, 1/4000,
+ 9, 7, 9, 8, 8, 7, 8, 7,
+ 1/9, 1/6, 1/8, 1/5), nrow = 4)
> DD <- D^2
> m <- apply(DD, 2, sum)
> m <- sqrt(m)
> M <- matrix(c(m, m, m, m), nrow = 4, byrow = T)
> R <- D/M
> R
[,1] [,2] [,3] [,4]
[1,] 0.3806169 0.5427204 0.5321521 0.3590803
[2,] 0.6343615 0.4221159 0.4656331 0.5386205
[3,] 0.4757711 0.5427204 0.5321521 0.4039654
[4,] 0.4757711 0.4824182 0.4656331 0.6463446
下面用資訊熵計算權重
> xxs <- function(x) -sum(x*log(x))/log(4)
> E <- apply(R, 2, xxs)
> F <- 1 - E
> w <- F/sum(F)
> w
[1] 0.32807915 0.10223509 0.04381521 0.52587055
于是權重為 0.32807915 0.10223509 0.04381521 0.52587055.
下面用Topsis法評價
> # 理想數法
> W <- diag(w)
> V <- R %*% W
> v_max <- apply(V, 2, max)
> v_min <- apply(V, 2, min)
> V_MAX <- matrix(c(v_max, v_max, v_max, v_max), nr = 4, byrow = TRUE)
> V_MIN <- matrix(c(v_min, v_min, v_min, v_min), nr = 4, byrow = TRUE)
> fun <- function(x) sqrt(sum(x^2))
> s_max <- apply(V-V_MAX, 1, fun)
> s_min <- apply(V-V_MIN, 1, fun)
> C <- s_min/(s_max + s_min)
> C
[1] 0.06842871 0.68438747 0.23006165 0.74631797
于是各目的地得分為0.06842871 0.68438747 0.23006165 0.74631797,因此選A4最佳.
- 檔案《民航里程數》給出了我國的民航里程資料(見附錄 1,民航里程資料),已知民航里程數Y與時間t之間的關系可以近似使用關系式 \(y = k - \frac{k^2}{(\frac{a}{x})^b-2}\)刻畫 , 其中k(5<k<100),a(0<a<1),b(0<b<1)為待估引數.請根據樣本使用最小二乘法估計這三個引數.
解:利用R求解
> # -- 第四題--
> library(pracma)
> x <- scan("clipboard")
Read 20 items
> f4 <- function(u){
+ f <- numeric(20)
+ for (i in 1:20){
+ f[i] <- x[i] - (u[1]-u[1]^2/((u[2]/i)^u[3]-2))
+ }
+ f
+ }
> sol4 <- lsqnonlin(f4, x0 = c(10, 1, 1))
> sol4
$x
[1] 20.0 0.8 0.6
$ssq
[1] 7.087413e-13
$ng
[1] 1.371535e-05
$nh
[1] 4.710978e-08
$mu
[1] 0.7030284
$neval
[1] 23
$errno
[1] 1
$errmess
[1] "Stopped by small x-step."
于是 \(a = 0.8,\ b = 0.6, \ k=20\).
- 某地區進行了用電量調查,得到的資料(見附錄 2,高峰用電量與總用電量資料),其中X代表每月用電量,Y 代表高峰時期每小時的用電量,
(1)利用樣本資料, 以 Y 為因變數,X 為自變數建立線性回歸方程擬合兩者的關系,并對所得
回歸方程進行顯著性檢驗,
(2)若根據歷史資料估計, 7 月份的平均用電量為 2000,試估計該月份的用電高峰時期每小時用電量,并求其 95%置信區間.
解:
(1)利用R求解
> # -- 第五題 --
> t5 <- read.table("clipboard", header = T)
> t5
X Y
1 679 79
2 292 44
3 1012 56
4 493 79
5 582 270
6 1156 364
7 997 473
8 2189 950
9 1097 534
10 2078 685
11 1818 584
12 1700 521
13 747 325
14 2030 443
15 1643 316
16 414 50
17 354 17
18 1276 188
19 745 77
20 435 139
21 540 56
22 874 156
23 1543 528
24 1029 64
25 710 400
26 1434 31
27 837 420
28 1748 488
29 1381 348
30 1428 758
31 1255 263
32 1777 499
33 370 59
34 2316 819
35 1130 479
36 463 51
37 770 174
38 724 410
39 808 394
40 790 96
41 783 329
42 406 44
43 1242 324
44 658 214
45 1746 571
46 468 64
47 1114 190
48 413 51
49 1787 833
50 3560 1494
51 1495 511
52 2221 385
53 1526 393
> x <- t5$X
> y <- t5$Y
> lm(y~x) -> sol5
> summary(sol5)
Call:
lm(formula = y ~ x)
Residuals:
Min 1Q Median 3Q Max
-413.99 -82.75 -19.34 123.76 315.22
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -83.13037 44.16121 -1.882 0.0655 .
x 0.36828 0.03339 11.030 4.11e-15 ***
---
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 157.7 on 51 degrees of freedom
Multiple R-squared: 0.7046, Adjusted R-squared: 0.6988
F-statistic: 121.7 on 1 and 51 DF, p-value: 4.106e-15
#常數項未通過檢驗,置為0
> lm(y~0+x) -> sol51
> summary(sol51)
Call:
lm(formula = y ~ 0 + x)
Residuals:
Min 1Q Median 3Q Max
-418.58 -113.30 -58.11 87.54 377.89
Coefficients:
Estimate Std. Error t value Pr(>|t|)
x 0.31351 0.01678 18.69 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 161.5 on 52 degrees of freedom
Multiple R-squared: 0.8704, Adjusted R-squared: 0.8679
F-statistic: 349.2 on 1 and 52 DF, p-value: < 2.2e-16
于是
\[y = 0.31351x \]且檢驗統計量 \(F = 349.2, \ p << 0.01\),因此認為兩者具有顯著的線性關系,\(R^2 = 0.8704\) 說明擬合程度較好.
來看看擬合效果
> plot(x,y)
> abline(sol51)
(2) R 程式如下
> x0 <- data.frame(x = 2000)
> Pred <- predict(sol51, x0, interval = "confidence", level = 0.95)
> Pred
fit lwr upr
1 627.0274 559.7 694.3549
于是得到高峰用電量為627.0274,其置信區間為 (559.7, 694.3549).
- 求解下面二次規劃問題:
請寫出(lingo 或者 R 或者 matlab)程式和計算結果,
解:規劃問題用Lingo求解比較方便,考慮使用Lingo求解
model:
min = x1^2 + x2^2 - 2*x1 - 5*x2;
x1 - 2*x2 >= -2;
-x1 - 2*x2 >= -6;
-x1 + 2*x2 >= -2;
end
輸出:
Global optimal solution found.
Objective value: -6.450000
Objective bound: -6.450000
Infeasibilities: 0.000000
Extended solver steps: 0
Total solver iterations: 8
Elapsed runtime seconds: 0.11
Model is convex quadratic
Model Class: QP
Total variables: 2
Nonlinear variables: 2
Integer variables: 0
Total constraints: 4
Nonlinear constraints: 1
Total nonzeros: 8
Nonlinear nonzeros: 2
Variable Value Reduced Cost
X1 1.400002 0.000000
X2 1.700001 0.000000
Row Slack or Surplus Dual Price
1 -6.450000 -1.000000
2 0.000000 -0.8000000
3 1.199997 0.000000
4 4.000000 0.000000
于是目標最小值為-6.45,其中 \(x_1 = 1.4,x_2=1.7\).
- 考慮 3 個產地與 4 個銷地的問題,各產地的產量、各銷地的銷量及各產地到各銷地的單 位運費如下表,問如何調運,才可使總的運費最少?請根據問題建立數學模型,
解:7. 發現產銷不平衡,產量(76)>銷量(66),因此假想一個銷地 B5,各產地到B5的運費為0且銷量為10.
設Ai到Bj的運量為 \(x_{ij}\),運費矩陣為
\[A = \left[ \begin{array}{c} 6&2&6&7&0\\ 4&9&5&3&0\\ 8&8&1&5&0 \end{array} \right] \]使用Lingo求解
代碼:
model:
sets:
thr /1..3/ : cl;
fiv /1..5/ : xl;
coor(thr,fiv) : A, x;
endsets
data:
cl = 30, 25, 21;
xl = 15, 17, 22, 12, 10;
A = 6, 2, 6, 7, 0
4, 9, 5, 3, 0
8, 8, 1, 5, 0;
enddata
min = @sum(coor(i, j): A(i,j)*x(i,j));
@for(thr(i) : @sum(fiv(j) : x(i,j)) = cl(i));
@for(fiv(j) : @sum(thr(i) : x(i,j)) = xl(j));
end
輸出:
Global optimal solution found.
Objective value: 161.0000
Infeasibilities: 0.000000
Total solver iterations: 8
Elapsed runtime seconds: 0.05
Model Class: LP
Total variables: 15
Nonlinear variables: 0
Integer variables: 0
Total constraints: 9
Nonlinear constraints: 0
Total nonzeros: 42
Nonlinear nonzeros: 0
Variable Value Reduced Cost
CL( 1) 30.00000 0.000000
CL( 2) 25.00000 0.000000
CL( 3) 21.00000 0.000000
XL( 1) 15.00000 0.000000
XL( 2) 17.00000 0.000000
XL( 3) 22.00000 0.000000
XL( 4) 12.00000 0.000000
XL( 5) 10.00000 0.000000
A( 1, 1) 6.000000 0.000000
A( 1, 2) 2.000000 0.000000
A( 1, 3) 6.000000 0.000000
A( 1, 4) 7.000000 0.000000
A( 1, 5) 0.000000 0.000000
A( 2, 1) 4.000000 0.000000
A( 2, 2) 9.000000 0.000000
A( 2, 3) 5.000000 0.000000
A( 2, 4) 3.000000 0.000000
A( 2, 5) 0.000000 0.000000
A( 3, 1) 8.000000 0.000000
A( 3, 2) 8.000000 0.000000
A( 3, 3) 1.000000 0.000000
A( 3, 4) 5.000000 0.000000
A( 3, 5) 0.000000 0.000000
X( 1, 1) 2.000000 0.000000
X( 1, 2) 17.00000 0.000000
X( 1, 3) 1.000000 0.000000
X( 1, 4) 0.000000 2.000000
X( 1, 5) 10.00000 0.000000
X( 2, 1) 13.00000 0.000000
X( 2, 2) 0.000000 9.000000
X( 2, 3) 0.000000 1.000000
X( 2, 4) 12.00000 0.000000
X( 2, 5) 0.000000 2.000000
X( 3, 1) 0.000000 7.000000
X( 3, 2) 0.000000 11.00000
X( 3, 3) 21.00000 0.000000
X( 3, 4) 0.000000 5.000000
X( 3, 5) 0.000000 5.000000
Row Slack or Surplus Dual Price
1 161.0000 -1.000000
2 0.000000 -2.000000
3 0.000000 0.000000
4 0.000000 3.000000
5 0.000000 -4.000000
6 0.000000 0.000000
7 0.000000 -4.000000
8 0.000000 -3.000000
9 0.000000 2.000000
于是方案為:
產地/銷地 | B1 | B2 | B3 | B4 | 產量 |
---|---|---|---|---|---|
A1 | 2 | 17 | 1 | 0 | 30 |
A2 | 13 | 0 | 0 | 12 | 25 |
A3 | 0 | 0 | 21 | 0 | 21 |
銷量 | 15 | 17 | 22 | 12 |
費用最小為161.
- 某工廠生產 A,B,C 三種產品,其所需勞動力、材料等有關資料如下表所示.(1)確定獲利最 大的生產方案;(2)產品 A,B,C 的利潤分別在什么范圍內變動時,上述最優方案不變;(3)如果勞動 力數量不增,材料不足可以從市場購買,每單位0.4元,問該廠要不要購進原材料擴大生產,購買 多少為宜?(要求給出數學模型以及相關的程式 lingo 或者 R 或者 matlab)
解:(1)設生產A,B,C產量分別為 \(x_1,x_2,x_3\) , \(z\) 為目標,可得如下模型:
\[\begin{aligned} &\max \quad z=3x_1+2x_2+4x_3\\ & \quad\quad6x_1+3x_2+5x_3 \le 45\\ & \quad\quad3x_1+4x_2+5x_3 \le 30\\ & \quad\quad x_1,x_2,x_3\ge0; \end{aligned} \]使用Lingo求解
代碼:
model:
max = 3*x1+2*x2+4*x3;
6*x1+3*x2+5*x3<=45;
3*x1+4*x2+5*x3<=30;
end
輸出結果為:
Global optimal solution found.
Objective value: 27.00000
Infeasibilities: 0.000000
Total solver iterations: 2
Elapsed runtime seconds: 0.05
Model Class: LP
Total variables: 3
Nonlinear variables: 0
Integer variables: 0
Total constraints: 6
Nonlinear constraints: 0
Total nonzeros: 12
Nonlinear nonzeros: 0
Variable Value Reduced Cost
X1 5.000000 0.000000
X2 0.000000 1.000000
X3 3.000000 0.000000
Row Slack or Surplus Dual Price
1 27.00000 1.000000
2 0.000000 0.2000000
3 0.000000 0.6000000
4 5.000000 0.000000
5 0.000000 0.000000
6 3.000000 0.000000
故最優方案為:生產5個A,0個B,3個C,且獲利最大為27.
(2) 靈敏度分析可得
Ranges in which the basis is unchanged:
Objective Coefficient Ranges:
Current Allowable Allowable
Variable Coefficient Increase Decrease
X1 3.000000 1.800000 0.6000000
X2 2.000000 1.000000 INFINITY
X3 4.000000 1.000000 1.000000
Righthand Side Ranges:
Current Allowable Allowable
Row RHS Increase Decrease
2 45.00000 15.00000 15.00000
3 30.00000 15.00000 7.500000
則A的范圍為[2.4,4.8],B的范圍為[0.0,3.0],C的范圍為[3.0,5.0],最優方案不變.
(3) 設購買原材料數量為\(x_4\).
模型改為
使用Lingo求解
代碼:
model:
max = 3*x1+2*x2+4*x3-0.4*x4;
6*x1+3*x2+5*x3 <= 45;
30+x4 = 3*x1+4*x2+5*x3;
end
輸出結果為:
Global optimal solution found.
Objective value: 30.00000
Infeasibilities: 0.000000
Total solver iterations: 1
Elapsed runtime seconds: 0.04
Model Class: LP
Total variables: 4
Nonlinear variables: 0
Integer variables: 0
Total constraints: 7
Nonlinear constraints: 0
Total nonzeros: 15
Nonlinear nonzeros: 0
Variable Value Reduced Cost
X1 0.000000 0.6000000
X2 0.000000 0.8000000
X3 9.000000 0.000000
X4 15.00000 0.000000
Row Slack or Surplus Dual Price
1 30.00000 1.000000
2 0.000000 0.4000000
3 0.000000 -0.4000000
4 0.000000 0.000000
5 0.000000 0.000000
6 9.000000 0.000000
7 15.00000 0.000000
故購買原材料15個最合適,利潤最大為30.
轉載請註明出處,本文鏈接:https://www.uj5u.com/houduan/485738.html
標籤:R