徐泽水《不确定多属性决策方法与应用》81页

第一步: 把原始决策矩阵A 利用适当的方法进行规范化为R，R为归一化后的矩阵

#### 第一步: 把原始决策矩阵A 利用适当的方法进行规范化为R，R为归一化后的矩阵
### norm_matrix()函数，根据书中收益型属性（按公式1.2）与成本型属性(按公式1.4)分别进行归一化
#####  注意这个与前面的norm_matrix函数(即下面norm_matrix2)结果相同，只是代码显得更少了,保证了列名不变
norm_matrix = function(A, shouyi = NULL, chengben = NULL) {
  stopifnot(!is.null(shouyi) | !is.null(chengben))
  if (is.matrix(A)) A = data.table(A)
  m = ncol(A)
  if (is.null(chengben)) chengben = setdiff(1:m, shouyi)
  if (is.null(shouyi)) shouyi = setdiff(1:m, chengben)
  # 如果输入的shouyi与chengben向量交集不为空，且并集不是全集，则算法出错
  stopifnot(length(intersect(shouyi, chengben)) == 0, setequal(union(shouyi, chengben), 1:m))
  R =copy(A) # 重新赋值
  if (length(chengben) == 0) {
    R[, colnames(R)[shouyi] := lapply(.SD, function(x) x / max(x)), .SDcols = shouyi] # 收益型属性归一化 （书中1.2式）
  } else if (length(shouyi) == 0) {
    R[, colnames(R)[chengben] := lapply(.SD, function(x) min(x) / x  ), .SDcols = chengben]# 成本型属性归一化 （书中1.3式）
  } else{
    R[, colnames(R)[shouyi] := lapply(.SD, function(x) x / max(x)), .SDcols = shouyi] # 收益型属性归一化 （书中1.2式）
    R[, colnames(R)[chengben] := lapply(.SD, function(x) min(x) / x  ), .SDcols = chengben]# 成本型属性归一化 （书中1.3式）
  }
  return(setDF(R))
}

norm_matrix2 = function(A,shouyi=NULL,chengben=NULL){
  if(is.matrix(A))A = data.table(A)
  stopifnot(!is.null(shouyi) | !is.null(chengben))
  m = ncol(A)
  if(is.null(chengben)) chengben =setdiff(1:m,shouyi) 
  if(is.null(shouyi)) shouyi = setdiff(1:m,chengben)
  stopifnot(length(intersect(shouyi,chengben))==0,setequal(union(shouyi,chengben),1:m))
  
  if( length(chengben) == 0 ){
    # 对决策矩阵进行重命名
    colnames(A)=paste0('V',1:m)
    shouyi = paste0("V",shouyi)
    R = A    
    # 归一化
    R[,':='(c(shouyi),lapply(.SD, function(x)x/max(x))),.SDcols =shouyi] # 收益型属性归一化 （书中1.2式）
  }else if( length(shouyi) == 0 ){
    #对决策矩阵进行重命名
    names(A)=paste0('V',1:m)
    chengben = paste0("V",chengben)
    R = A
    # 归一化
    R[,':='(c(chengben),lapply(.SD,function(x)min(x)/x)),.SDcol = chengben]# 成本型属性归一化 （书中1.3式）
  }else{
    #对决策矩阵进行重命名
    names(A)=paste0('V',1:m)
    shouyi = paste0("V",shouyi)
    chengben = paste0("V",chengben)
    R = A
    # 归一化
    R[,':='(c(shouyi),lapply(.SD, function(x)x/max(x))),.SDcols =shouyi] # 收益型属性归一化 （书中1.2式）
    R[,':='(c(chengben),lapply(.SD,function(x)min(x)/x)),.SDcol = chengben]# 成本型属性归一化 （书中1.3式）    
  }
  R = as.data.frame(R)
  return(R)
}

R = norm_matrix(A, chengben = 4)
round(R,3)
#>       V1    V2    V3    V4    V5
#> 1  0.799 1.000 1.000 0.723 1.000
#> 2  0.734 0.547 0.411 0.754 0.535
#> 3  1.000 0.833 0.682 0.846 0.816
#> 4  0.793 0.638 0.395 1.000 0.470
#> 5  0.706 0.797 0.522 0.923 0.592
#> 6  0.448 0.612 0.268 0.838 0.625
#> 7  0.650 0.721 0.539 0.849 0.675
#> 8  0.979 0.620 0.511 0.977 0.585
#> 9  0.489 0.462 0.238 0.723 0.574
#> 10 0.658 0.537 0.380 0.875 0.554
#> 11 0.520 0.654 0.467 0.740 0.725
#> 12 0.421 0.648 0.272 0.731 0.721
#> 13 0.456 0.562 0.344 0.746 0.687
#> 14 0.394 0.497 0.290 0.690 0.547
#> 15 0.479 0.489 0.357 0.768 0.581
#> 16 0.366 0.430 0.524 0.635 0.715

第二步 ： 求出综合属性正理想值z_max,以及综合属性负理想值z_min

###### 第二步 ： 求出综合属性正理想值z_max,以及综合属性负理想值z_min
M_zonghe_position = function(R,lower_c,upper_c){
   #综合属性正理想值z_max
    library(Rglpk)
    n = nrow(R)
    m = ncol(R)
    ## 约束条件，权和向量为1 
    mat = matrix(rep(1,m),nrow = 1) 
    dir = c("==")
    rhs = c(1)
    ## 
    types = c("C") # 表示解为实数
    bounds <- list(lower = list(ind = 1L:m, val = lower_c),
                   upper = list(ind = 1L:m, val = upper_c))
    ###  下面 max_obj 函数中的xx为R矩阵的某一行
    max_obj = function(xx) Rglpk_solve_LP(xx, mat, dir, rhs, bounds, types,max = TRUE)$optimum
    return(apply(R, 1, max_obj))
}
z_max = M_zonghe_position(R,lower_c = c(0.22,0.18,0.15,0.23,0.16),upper_c = c(0.24,0.20,0.17,0.26,0.17))
round(z_max,3)
#>  [1] 0.890 0.623 0.851 0.706 0.735 0.585 0.703 0.777 0.522 0.633 0.631 0.576
#> [13] 0.575 0.502 0.555 0.534


## 求出综合属性负理想值
M_zonghe_negative = function(R,lower_c,upper_c){
  #综合属性负理想值z_min
  if(is.data.table(R)) R = as.data.frame(R)
  library(Rglpk)
  n = nrow(R)
  m = ncol(R)
  ## 约束条件，权和向量为1 
  mat = matrix(rep(1,m),nrow = 1) 
  dir = c("==")
  rhs = c(1)
  ## 
  types = c("C") # 表示解为实数
  bounds <- list(lower = list(ind = 1L:m, val = lower_c),
                 upper = list(ind = 1L:m, val = upper_c))
  ###  下面 max_obj 函数中的xx为R矩阵的某一行
  min_obj = function(xx) Rglpk_solve_LP(xx, mat, dir, rhs, bounds, types)$optimum
  return(apply(R, 1, min_obj))
}
z_min = M_zonghe_negative(R,lower_c = c(0.22,0.18,0.15,0.23,0.16),upper_c = c(0.24,0.20,0.17,0.26,0.17))
z_min
#>  [1] 0.8799025 0.6122816 0.8442960 0.6869648 0.7224206 0.5680731 0.6945892
#>  [8] 0.7600534 0.5084315 0.6185512 0.6213819 0.5605654 0.5620055 0.4893540
#> [15] 0.5430971 0.5233214

第三步： 求出各方按的满意度，PW 矩阵的每一行为对应方案的满意度矩阵

#####  第三步： 求出各方按的满意度，PW 矩阵的每一行为对应方案的满意度矩阵，
# 把满意度矩阵的每一行乘以对应属性的权重，则为该方案的满意度，但此时属性权重未知
# 通过建立线性目标函数，求解出属性权重
PW = matrix(0,ncol = ncol(R),nrow = nrow(R))
for(i in 1:nrow(R)){
  for(j in 1:ncol(R)){
    PW[i,j] = (R[i,j] - z_min[i])/(z_max[i] - z_min[i])
  }
}
PW
#>             [,1]        [,2]         [,3]        [,4]       [,5]
#>  [1,]  -7.820666 11.65256149  11.65256149 -15.1896961  11.652561
#>  [2,]  11.145778 -6.00900673 -18.47012741  12.9486991  -7.058863
#>  [3,]  22.923509 -1.62874601 -23.94335797   0.2436666  -4.150187
#>  [4,]   5.606157 -2.60510136 -15.40681623  16.5088261 -11.448630
#>  [5,]  -1.373481  6.09491123 -16.36028397  16.3793613 -10.672955
#>  [6,]  -7.029567  2.55423994 -17.59326739  15.8143586   3.311256
#>  [7,]  -4.979815  2.92839875 -17.50474858  17.3674915  -2.235035
#>  [8,]  13.013027 -8.33758206 -14.79304717  12.8414038 -10.367462
#>  [9,]  -1.441419 -3.43690273 -20.14420736  15.9848337   4.883279
#> [10,]   2.696916 -5.63428706 -16.47249835  17.6965044  -4.490027
#> [11,] -10.399110  3.40588154 -15.93163587  12.2046922  10.724431
#> [12,]  -9.130713  5.74617525 -18.86296457  11.1727327  10.494447
#> [13,]  -7.990561  0.02312869 -16.44500376  13.8874395   9.466552
#> [14,]  -7.630903  0.58843249 -15.97617926  16.0642115   4.593201
#> [15,]  -5.286466 -4.41756470 -15.36276214  18.5132509   3.161791
#> [16,] -15.310838 -9.07672150   0.08382542  10.8811429  18.581452

### 通过满意度矩阵，求出属性权重
M_w = function(R,lower_c,upper_c){
  library(Rglpk)
  n = nrow(R)
  m = ncol(R)
  obj = c(apply(R, 2, sum)) # 设置目标函数
  mat = matrix(rep(1,m),nrow = 1) # 约束条件，权和向量为1 
  dir = c("==")
  rhs = c(1)
  types = c("C")
  bounds <- list(lower = list(ind = 1L:m, val = lower_c),
                 upper = list(ind = 1L:m, val = upper_c))
  return(Rglpk_solve_LP(obj, mat, dir, rhs, bounds, types,max =TRUE )$solution)
}
w = M_w(PW,lower_c = c(0.22,0.18,0.15,0.23,0.16),upper_c = c(0.24,0.20,0.17,0.26,0.17))
w # 权重
#> [1] 0.22 0.20 0.15 0.26 0.17
### 通过权重确定方案综合指标值
z = apply(R, 1, function(x)sum(x*w))
z
#>  [1] 0.8839165 0.6193408 0.8475817 0.7011357 0.7335354 0.5851453 0.7034886
#>  [8] 0.7693701 0.5212981 0.6306097 0.6310837 0.5758520 0.5752496 0.5018189
#> [15] 0.5552379 0.5317315

##第四步，#按降序排列，最大的为方案最优

#########第四步，#按降序排列，最大的为方案最优

round(z,4)%>% rank %>% order(.,decreasing=T)#按降序排列，最大的为最优
#>  [1]  1  3  8  5  7  4 11 10  2  6 12 13 15 16  9 14