求教一个优化问题

x5miao

别人问我的一个问题，开始觉得很简单，但考虑了很长时间发现其计算过程还是很有难度的。

题目: 有三相电A,B,C，分别有n,m,k个用户（电表），各个相限内各个住户在一天24小时内的用电量不尽相同，造成三相电用电不平衡。电力部门希望对辖区内所有相限内住户进行调整，使调整之后的各个相限用电量差值最小。EAit表示A相限第i个用户在t时间内的用电量，以此类推。i=1.2.3,……n,t=1,2,3……24.要求用数学符号进行推导和表达，均可假设条件

我自己分析了一下，分析过程如下：
假设1：用户在每个时刻使用的电量与其所使用的电的象限无关（不然题意不足无法求解）；
假设2：三相电的供电能力是相等的，不对用户的使用产生限制； 假设3：我们使用ABC三相电在每个时刻输出电量的差值的平方和在一天之内的累加值大小作为衡量“各个象限用电量差值”的依据。也即


接下来怎么求就不知道了，请问有没有比穷举更高效的算法来解决这个问题？


A[5][24] = {
{0,0,0,0,0,3,8,5,4,6,10,14,12,5,5,7,4,10,16,14,13,9,7,1},
{ 3,3,4,2,0,0,2,1,1,0,2,9,2,1,1,9,16,12,16,8,12,2,4,2},
{10,11,14,12,10,10,5,2,9,14,20,8,2,10,10,18,,10,10,10,2,10,10,10,17},
{4,3,6,3,1,2,9,11,2,4,1,9,8,7,10,15,2,4,,6,2,3,5,2,1,0},
{10,11,14,0,4,10,5,2,1,3,7,0,0,12,3,18,5,8,5,7,10,12,0,0},
};
B[6][24] = {
{10,11,4,12,10,30,5,2,9,14,20,8,2,10,10,18,10,1,10,2,11,10,2,7},
{0,11,14,12,4,10,5,2 ,9,14,20,8,2,10,10,10,10,10,10,2,10,3,2,1},
{0,1,14,14,10,4,5,2,9,14,20,8,2,10,10,18,10,10,10,2,10,10,10,17},
{0,7,14,12,3,10,5,0,9,0,0,0,2,2,10,5,10,4,4,2,1,0,0,0},
{2,3,14,12,10,10,5,2,9,14,20,8,2,10,10,18,10,10,10,2,10,10,10,17},
{28,11,14,12,10,10,5,2,9,14,20,8,2,10,10,18,10,10,10,2,10,10,10,17},
};
C[4][24] = {
{0,11,0,0,0,0,5,2,12,1,13,4,2,10,2,18,10,10,10,2,10,10,10,17},
{0,0,14,12,10,10,5,2,9,14,20,2,2,2,10,18,10,10,10,2,0,0,0,0},
{0,11,0,12,4,4,5,2,9,14,20,12,2,11,10,18,10,10,10,2,10,2,1,1},
{10,4,14,12,10,10,5,2,9,14,20,8,2,10,10,18,10,10,10,2,10,10,0,0},
};

tangboyun

你的objFunc写复杂了，因为每个小时的总输出功率是事先就知道的，所以我觉得不需要互相减，每个组之间相互接近，也就是都接近每小时总功率的1/3。求这个mse即可，这样简单很多。如此一来objFunc就不需要考虑不同维度的矩阵运算了（col为24,row为n、m、k）。接下来由于每个用户各带自己那24维的功率向量，所以可以把每个用户和各自带的那个24向量和在一起作为一个变量来考虑，因为最终objFunc里不需要你那样互相减了，直接一个矩阵减法接平方和就是了。
问题简化后就变成有一共有n +m +k个不同用户，最后需要分成三组个数分别为m、n、k，然后min objFunc

我个人觉得可以用拉格郎日乘子法求偏导后很容易的解决掉。

上面这句错了。我明天再想想。

fallening

这不是一个整数规划么？

x5miao

回复 2# tangboyun


    o，那个M在某一固定时刻是常量，却是对时间的变量。

tangboyun

to楼上，我知道啊，但是你想想objFunc里，只要考虑24个小时这24个离散的量，每个用户每小时各消耗多少都清楚的啊，mse的话，只要求 1 到 24小的每组减去该小时总消耗量（已经知道）/3的mse之和啊。objFunc里，只要划组定了。最后的结果就定了。
minimize      
       for 1 to 24  { sum (Atotal - 1/3该小时总量)^2 + (Btotal - 1/3 该小时总量)^2 + (Ctotal - 1/3 该小时总量)^2 }

不需要相互减了把，这样简单很多。

tangboyun

初步设想如下：
假设我5楼的那个objFunc正确的话，那么可以以如下基于swap的方式求得最优解。
每个用户为1x24向量，总用户的输出功率为（n+m+k）×24的矩阵，则每组平均 Garvg为将上述矩阵列求和后每个元素除以组数的1x24向量
最优值objFunc为求各个组的Sum Square Error的最小值。(Atotal - Garvg)^2 + (Btotal - Garvg)^2 + (Btotal-Garvg)^2

初始化时，将用户按个数随机划入各组。以下为递归迭代步骤：

1、按各组的SQE，将各组降序排序，令排序结果为G1,G2,G3。
2、算法主要不停的迭代选取1对最有可能改善当前状况的用户进行组间交换(令为c1、c2)。选取的顺序保证从G1挑选第一个用户，然后依次尝试G2、G3（挑第二个）。如果选取成功，更新组间信息，回到步骤1继续迭代。直到各组尝试皆失败，算法终止。（我不知道需不需要尝试G2、G3这种组合，没数学证明，如果都尝试肯定能拿到最优解）
      挑选的标准如下在选中的2个组中选取能使该值最小的一对（G1total - Garvg - （c1 -c2））^2 + (Gxtotal - Garvg + (c1-c2))^2
     关于（c1-c2）这个可以事先建（m+n+k) × (m+n+k)的上三角矩阵每个角标存这个向量。

关于算法稳定性的测试：
       可以以相同数据不同随机数种子划分比如100次实验，如果最后的最优值相同。则说明能逼近最优解。需要说明的是，最优方案肯定不只一种（比如1 2 1的这种划分，两个1用户的组，最后可以互换），但最优值应该只有一个。

cjaizss

此类背包,除了贪婪求近似解之外,你不穷举还想怎么办

tangboyun

用Haskell写了一个解这个问题的小程序，带测试

编译方法:
ghc --make -O3 test.hs
运行：
3、4、5指分别生成3个组，每个用户为3、4、5，该程序可以指定任意组，任意用户数。当然你要把24改成其他也行。
./test 3 4 5

自动roll随机数初始化，取得数据后求解一次，然后随机分割后求解100次，并比较结果稳定性。

发现无法保证每次都得到最优解。只能求得一个局部最优。代码写的应该比较容易读。有兴趣自己尝试把。

1. 
   
2. -- file test.hs
   
3. {-# LANGUAGE TupleSections #-}
   
4. module Main where
   
5. 
   
6. import qualified Data.Vector.Unboxed as UV
   
7. import qualified Data.Vector as V
   
8. import Data.List
   
9. import Data.Ord
   
10. import Data.Array
    
11. import System.Random
    
12. import Control.Monad
    
13. import Data.Time.Clock
    
14. import System.Environment
    
15. 
    
16. data SplitTable = S {
    
17.   groupId :: String,
    
18.   groupMember :: [Int],
    
19.   gapFromGroupSQE :: Double,  -- ^ 与均值向量的sum square error,排序准则
    
20.   gapVec :: Data              -- ^ Gtotal - Gmean
    
21.   } deriving (Eq)
    
22. 
    
23. -- For Sort
    
24. instance Ord SplitTable where
    
25.   compare = comparing gapFromGroupSQE
    
26. -- For Pretty Print  
    
27. instance Show SplitTable where
    
28.   show a = (groupId a) ++ ": " ++
    
29.            (unwords $ map show $ groupMember a)
    
30.            
    
31. type Users = V.Vector Data
    
32. 
    
33. type Data = UV.Vector Double
    
34. type GapTable = Array (Int,Int) (UV.Vector Double)
    
35. 
    
36. atU :: UV.Unbox a => UV.Vector a -> Int -> a
    
37. atU = UV.unsafeIndex
    
38. at :: V.Vector a -> Int -> a
    
39. at = V.unsafeIndex
    
40. 
    
41. calcSQE :: Users -> Data -> Double
    
42. calcSQE usrs meanVec =
    
43.   UV.sum $ UV.map (^2) $ (V.foldl1' (.+) usrs) .- meanVec
    
44.   
    
45. calcGap :: Users -> Data -> Data  
    
46. calcGap usrs meanVec =
    
47.   V.foldl1' (.+) usrs .- meanVec
    
48. 
    
49. -- pairwise add
    
50. (.+) :: Data -> Data -> Data
    
51. (.+) = UV.zipWith (+)
    
52. -- pairwise minus
    
53. (.-) :: Data -> Data -> Data
    
54. (.-) = UV.zipWith (-)
    
55. -- pairwise divide  
    
56. (./.) :: Data -> Double -> Data
    
57. vec ./. d = UV.map (/ d) vec      
    
58. 
    
59. initialization :: StdGen    -- ^ 随机数生成器
    
60.                -> Int    -- ^ 小时数，设短的方便测试
    
61.                -> [Int]  -- ^ 各组的用户数
    
62.                -> (Users,GapTable)
    
63. initialization gen hourPerDay numList =
    
64.   let groupVec = UV.fromList numList
    
65.       groupNum = UV.length groupVec
    
66.       totalUser = UV.sum groupVec
    
67.       rndNumbers = UV.fromList $ take (totalUser * hourPerDay) $
    
68.                    randomRs
    
69.                    (0.0,20.0) -- ^ 功率的范围
    
70.                    gen
    
71.       allUsers = V.fromList $
    
72.                  map (\idx ->
    
73.                        UV.slice (hourPerDay * idx) hourPerDay rndNumbers
    
74.                      ) [0..totalUser - 1]
    
75.       gTable    = array ((0,0),(totalUser-1,totalUser-1))
    
76.                   [((idx1,idx2),(allUsers `at` idx1) .-
    
77.                                 (allUsers `at` idx2))
    
78.                   |idx1 <- [0..totalUser-1],idx2 <- [0..totalUser-1]]
    
79.   in (allUsers,gTable)
    
80. 
    
81.    
    
82. randomSplitUsers :: StdGen  -- ^ 随机数生成器
    
83.                  -> [Int]
    
84.                  -> Users
    
85.                  -> [SplitTable]
    
86. randomSplitUsers gen numList usrs =   
    
87.   let
    
88.     len = V.length usrs
    
89.     rs = UV.fromList $ take len $
    
90.          randoms $ gen :: UV.Vector Int
    
91.     shuffledIdx = UV.fromList $ sortBy (comparing (rs `atU`)) [0..len-1]
    
92.     groupVec = UV.fromList numList
    
93.     groupMean = (V.foldl1' (.+) usrs) ./. fromIntegral groupNum
    
94.     groupNum = UV.length groupVec
    
95.     slVec = UV.zip (UV.prescanl' (+) 0 groupVec) groupVec
    
96.   in map (\idx ->
    
97.            let (begIdx,len) = slVec `atU` idx
    
98.                ls = sort $ UV.toList $
    
99.                     UV.slice begIdx len shuffledIdx
    
100.                selectd = V.ifilter (\i _ ->
     
101.                                      i `elem` ls) usrs
     
102.            in
     
103.             S {
     
104.               groupId = "Group" ++ show (idx+1),
     
105.               groupMember = ls,
     
106.               gapFromGroupSQE = calcSQE selectd groupMean,
     
107.               gapVec = calcGap selectd groupMean
     
108.               }
     
109.          ) [0..groupNum-1]
     
110.    
     
111. -- | 目标为最小化该函数
     
112. objFunc :: [SplitTable] -> Double
     
113. objFunc = sum . map gapFromGroupSQE
     
114. 
     
115. 
     
116. 
     
117. evolve :: Users -> GapTable -> [SplitTable] -> [SplitTable]
     
118. evolve allUsers gt ts = go (objFunc ts) (sort' ts)
     
119.   where
     
120.     groupMean = (V.foldl1' (.+) allUsers) ./. fromIntegral (length ts)
     
121.     sort' = sortBy (\a b -> compare b a) -- ^ 降序
     
122.     go objV ts' =
     
123.       let combs = concatMap (\ss -> map (head ss,) (tail ss)) $
     
124.                   filter ((>1).length) (tails ts')
     
125.           pair =  msum $ map selectPairIn combs -- ^ Monad Just搜到就终止
     
126.       in case pair of
     
127.         Just ((t1,i),(t2,j)) ->
     
128.           let
     
129.             (t1',t2') = swap ((t1,i),(t2,j))
     
130.             newTs = sort' $ t1':t2':(ts' \\ [t1,t2])
     
131.             newObjV = objFunc newTs
     
132.           in go newObjV newTs
     
133.         Nothing    -> ts'
     
134.       where
     
135.         selectPairIn :: (SplitTable,SplitTable)
     
136.                      -> Maybe ((SplitTable,Int),(SplitTable,Int))
     
137.         selectPairIn (table1,table2) =
     
138.           let
     
139.             idxs1 = groupMember table1
     
140.             idxs2 = groupMember table2
     
141.             gap1 = gapFromGroupSQE table1
     
142.             gap2 = gapFromGroupSQE table2
     
143.             gapV1 = gapVec table1
     
144.             gapV2 = gapVec table2
     
145.             f = UV.sum . UV.map (^2)
     
146.             ls = [((f $ gapV1 .- (gt ! (idx1,idx2))) +
     
147.                    (f $ gapV2 .+ (gt ! (idx1,idx2))) -
     
148.                    gap1 - gap2,(idx1,idx2))
     
149.                  |idx1<-idxs1,idx2<-idxs2,
     
150.                   let newG1 =  f $ gapV1 .- (gt ! (idx1,idx2))
     
151.                       newG2 =  f $ gapV2 .+ (gt ! (idx1,idx2))
     
152.                       gap = newG1 + newG2 - gap1 - gap2
     
153.                   in gap < 0]
     
154.           in if null ls
     
155.              then Nothing
     
156.              else
     
157.                let (i1,i2) = snd $ minimumBy (comparing fst) ls
     
158.                in Just ((table1,i1),(table2,i2))
     
159.         swap :: ((SplitTable,Int),(SplitTable,Int))
     
160.              -> (SplitTable,SplitTable)
     
161.         swap ((table1,idx1),(table2,idx2)) =
     
162.           let mem1 = idx2 : (groupMember table1 \\ [idx1])
     
163.               mem2 = idx1 : (groupMember table2 \\ [idx2])
     
164.               gapV1 = gapVec table1
     
165.               gapV2 = gapVec table2
     
166.               newGap1 = gapV1 .- (gt ! (idx1,idx2))
     
167.               newGap2 = gapV2 .+ (gt ! (idx1,idx2))
     
168.               newG1 =  UV.sum $ UV.map (^2) $ newGap1
     
169.               newG2 =  UV.sum $ UV.map (^2) $ newGap2
     
170.           in (table1 {
     
171.                  groupMember = mem1,
     
172.                  gapFromGroupSQE = newG1,
     
173.                  gapVec = newGap1
     
174.                  } ,
     
175.               table2 {
     
176.                  groupMember = mem2,
     
177.                  gapFromGroupSQE = newG2,
     
178.                  gapVec = newGap2
     
179.                  }
     
180.              )
     
181. 
     
182. 
     
183. 
     
184. main :: IO ()
     
185. main = do
     
186.   numList <- fmap (map read) getArgs
     
187.   str <- fmap (show.utctDayTime) getCurrentTime
     
188.   let seed = read $ takeWhile (/= '.') str
     
189.       (g1,g2) = split $ mkStdGen seed
     
190.       (users,gt) = initialization g1 24 numList      
     
191.       tables = randomSplitUsers g2 numList users
     
192.   putStrLn "User List:"     
     
193.   mapM_ (\idx ->
     
194.           -- round to 0.01
     
195.           let v = map ((/ 100.0).fromIntegral.round.(* 100)) $
     
196.                   UV.toList $ users `at` idx
     
197.           in
     
198.            putStrLn $ "User" ++ show idx ++ ": " ++ (unwords $
     
199.                    map show v)
     
200.           ) $ [0..(V.length users)-1]
     
201.   putStrLn "Split to: "  
     
202.   putStrLn $ "ObjValue: " ++ show (objFunc tables)
     
203.   mapM_ print tables
     
204.   let result = evolve users gt tables
     
205.       objV = objFunc result
     
206.   putStrLn "After optimization:"
     
207.   mapM_ print result
     
208.   putStrLn $ "ObjValue: " ++ show objV
     
209.   -- ^ test same data 100 times with different shuffle
     
210.   let  
     
211.     testResult = msum $
     
212.         map (\rndSeed ->
     
213.               let ts = randomSplitUsers (mkStdGen rndSeed) numList users
     
214.                   ts' = evolve users gt ts
     
215.                   newObjV = objFunc $ ts'
     
216.               in if abs(objV - newObjV) < 1e-6
     
217.                  then Nothing
     
218.                  else Just $ "Error:\nNewObjV Found: " ++ show newObjV ++
     
219.                       "\n" ++ (unlines $ map show ts')
     
220.               )  [1..100]
     
221.   case testResult of
     
222.     Nothing -> putStrLn "100 Tests passed."
     
223.     Just str -> putStrLn str
     

复制代码

tangboyun

下面是Roll 2次的结果。
[tangboyun@TangBoYun-NoteBook 12]$ ./Opt 3 2 4
User List:
User0: 7.33 15.24 8.25 14.13 10.94 13.32 7.42 7.94 4.61 1.08 1.89 10.74 1.63 8.57 5.29 1.49 13.98 1.16 5.07 0.59 19.47 3.0e-2 8.57 4.75
User1: 10.84 16.99 1.07 2.59 10.47 13.54 5.48 8.19 6.63 0.82 16.83 5.0e-2 19.35 15.41 3.09 12.06 11.29 7.33 11.68 1.59 1.39 0.77 14.26 16.41
User2: 4.01 10.12 3.32 14.64 8.24 5.37 18.14 19.47 12.37 2.43 14.52 9.59 1.45 3.81 2.31 4.57 9.43 4.51 9.72 3.5 18.08 17.72 3.5 19.53
User3: 15.02 15.89 15.24 2.0 4.38 3.0 15.66 17.15 8.92 8.0e-2 4.54 18.2 12.73 10.55 3.44 2.59 11.35 3.57 18.78 13.45 10.48 6.33 12.29 6.4
User4: 11.53 11.8 13.73 8.83 1.77 13.48 2.89 17.04 2.98 6.94 0.43 13.87 5.44 10.09 1.34 18.92 11.52 8.61 3.43 12.46 11.6 9.0e-2 2.76 3.28
User5: 8.0 9.05 9.54 8.59 7.99 2.1 19.21 2.87 0.67 13.57 0.74 6.9 0.18 10.91 1.58 11.01 13.87 13.1 2.05 9.52 5.32 18.4 10.37 7.56
User6: 15.68 8.84 18.05 9.22 16.63 3.72 12.74 9.4 15.39 8.18 19.69 3.58 14.88 4.66 3.95 13.16 13.7 2.87 17.6 17.77 11.77 2.07 11.64 10.75
User7: 1.9 2.44 4.23 8.13 9.2 2.41 9.12 8.57 4.76 1.92 1.29 5.66 9.14 9.02 2.67 0.36 5.9 10.16 9.77 7.73 6.87 10.63 13.74 14.98
User8: 2.0 19.61 19.85 0.4 8.93 18.86 2.64 18.64 19.36 19.91 12.93 0.49 1.99 7.6 4.03 4.51 11.12 9.38 8.8 10.07 9.36 18.91 19.49 2.92
Split to:
ObjValue: 8509.47738783689
Group1: 2 5 6
Group2: 3 7
Group3: 0 1 4 8
After optimization:
Group2: 8 3
Group3: 5 7 0 1
Group1: 4 2 6
ObjValue: 4981.239332351991
Error:
NewObjV Found: 5028.666369780052
Group3: 7 2 1 4
Group2: 8 3
Group1: 0 5 6

[tangboyun@TangBoYun-NoteBook 12]$ ./Opt 3 2 4
User List:
User0: 4.77 11.83 1.51 12.3 14.9 11.87 12.77 19.53 12.37 10.65 12.9 5.12 17.79 8.24 17.05 10.97 6.15 15.43 18.75 4.33 14.01 14.76 10.04 9.74
User1: 19.97 4.78 12.79 5.41 15.33 7.63 1.43 19.44 14.24 12.89 13.13 13.34 12.78 0.62 3.34 2.63 2.64 10.2 8.72 2.01 9.37 7.14 16.68 7.18
User2: 11.74 9.32 14.74 2.56 5.44 6.4 13.8 8.87 11.33 18.22 5.67 4.78 12.78 7.09 7.04 10.91 8.84 10.41 4.0e-2 16.63 7.53 13.1 19.63 0.36
User3: 1.02 17.43 17.99 2.04 12.71 7.97 11.95 7.78 15.18 17.07 4.11 16.01 4.51 13.48 2.77 10.03 0.72 18.15 15.12 10.17 5.0e-2 16.58 16.63 10.26
User4: 15.39 16.13 0.96 16.6 5.62 10.53 11.66 11.24 3.8 1.19 8.64 12.77 4.52 19.54 12.14 8.67 5.55 11.55 3.5 16.88 4.18 7.53 12.72 16.29
User5: 9.24 3.37 16.61 18.41 19.65 7.1 8.14 9.01 19.4 2.57 18.35 14.8 17.33 1.02 15.19 18.83 11.62 16.91 5.18 0.73 3.21 8.38 8.1 11.46
User6: 0.35 10.7 0.43 5.73 5.58 1.79 8.15 12.46 13.81 2.72 9.35 19.64 9.88 6.49 18.23 16.68 3.91 16.14 5.96 3.0e-2 18.77 17.64 2.16 1.33
User7: 2.68 3.84 17.54 9.82 16.36 17.27 15.53 12.51 8.02 18.17 5.52 5.87 12.16 12.68 11.56 9.54 3.57 10.61 1.04 9.48 14.55 15.74 6.47 7.85
User8: 15.05 3.23 8.96 3.07 2.73 2.75 2.31 4.92 18.09 10.68 18.66 18.71 2.23 15.52 15.06 6.77 15.81 10.94 14.73 1.97 1.54 13.14 17.29 18.99
Split to:
ObjValue: 10196.90324190425
Group1: 0 1 6
Group2: 2 3
Group3: 4 5 7 8
After optimization:
Group2: 5 3
Group3: 2 4 1 6
Group1: 7 8 0
ObjValue: 5725.679217106968
All Tests passed.