均匀设计表的构造+考虑子目标偏好度的均匀设计

   其中均匀设计表的构造权重向量，主要是均匀构造表的生成。其中网上有详细的构造均匀设计表的步骤。

用到均匀设计，是因为项目中（多目标进化问题）中设计子目标的权向量，以下是我在项目中写的代码。（C#）

解释：GetMvector（）是原本的均匀设计构造表的生成。GetMvectorOfP(double[] P,double delta)则是集成了子目标偏好度的均匀设计表的生成。

其中思想步骤如下：

根据用户对不同子目标的偏好<p1,p2,p3,p4>，p1>=p2>=p3>=p4，p1+p2+p3+p4=1,挑选n个子问题权向量

  

   1）产生N0=n*（1+delta1）个权向量，0<=delta1<=0.4，组成集合N(0)；

   2）对N0个权向量，分别按照各个维度上的权重值，从大到小进行排序，得到4种排序结果；

   3）依据第1个维度上的排序结果，从N(0)中挑选前n*p1个权向量，，组成集合N(1)，作为偏好子目标1的的权向量；

   4）依据第2个维度上的排序结果，从N(0)-N(1)中挑选前n*p2个权向量，组成集合N(2)，作为偏好子目标2的的权向量；

   5）类似地，按照第3、4维度上的排序结果，挑选偏好子目标3、4的权向量集合N(3)、N(4);

   6）N(1)、N(2)、N(3)、N(4)组成全部n个子问题的权向量集合。

 

源代码如下：

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace MOEA_PCI
{
    public class UniformDesign
    {
        public int N;//N行
        public int M;//M列
        public UniformDesign(int n, int m)
        {
            if (IsPrime(n))
            {
                N = n;
            }
            else
            {
                N = GetNearPrime(n);
            }
            M = m;
        }
        public bool IsPrime(int n)//判断一个数是否是素数
        {
            if (n == 1 || n == 2 || n == 3 || n == 5)
            {
                return true;
            }
            else
            {
                for (int i = 2; i < n; i++)
                {
                    if (n % i == 0)
                    {
                        return false;
                    }
                }
                return true;
            }
        }
        public bool CdivisorofXY(int x, int y)//判断两个数的最大公约数是否是1
        {
            for (int i = 2; i <= x && i <= y; i++)
            {
                if (x % i == 0 && y % i == 0)
                {
                    return false;
                }
            }
            return true;
        }
        public List<int> GetAllBaseofGeneratingvVector()//获得所有能够产生生成向量的基数,尽量N 是素数
        {
            List<int> r = new List<int>();//最后的基数列表
            List<int> hv = new List<int>();//全部与N最大公约数为1的数，从这些数中选择一些判断是否能产生生成向量A
            for (int i = 2; i < N; i++)
            {
                if (CdivisorofXY(i, N))
                {
                    hv.Add(i);
                }
            }
            for(int i=0;i<hv.Count;i++)
            {
                List<int>tempv=new List<int>();
                bool ishv = false;
                for (int j = 1; ;j++ )
                {
                    int temp=(int)(Math.Pow(hv[i], j) % N);
                    if (temp == 1 && j-1>=M-2)
                    {
                        ishv = true;
                        break;
                    }
                    if (!tempv.Contains(temp))
                    {
                        tempv.Add(temp);
                    }
                    else
                    {
                        break;
                    }
                }
                if (ishv)
                {
                    r.Add(hv[i]);
                }
            }
            return r;
        }
        public double[,] GetM_1vector()//获得所有生成的N*M-1矩阵中CD2值最小的一个矩阵
        {
            List<int> bv = GetAllBaseofGeneratingvVector();//所有基数列表
            List<List<int>> mhv = new List<List<int>>();//所有的生成向量列表
            int minindex = 0;
            double mincd2 = 0.0;
            for (int i = 0; i < bv.Count; i++)
            {
               List<int> t=new List<int>();
                for(int j=0;j<M-1;j++)
                {
                    int temp = (int)(Math.Pow(bv[i], j) % N);
                    t.Add(temp);
                }
                mhv.Add(t);
            }
            for (int bvindex = 0; bvindex < bv.Count; bvindex++)
            {
                double[,] vector;//N*M-1矩阵
                vector = new double[N, M - 1];//生成N*M-1 矩阵
                for (int i = 0; i < N; i++)//N
                {
                    for (int j = 0; j < M - 1; j++)//m-1
                    {
                        if (i == 0)
                        {
                            vector[i, j] = mhv[bvindex][j];
                        }
                        else
                        {
                            if (vector[i - 1, j] + mhv[bvindex][j] <= N)
                            {
                                vector[i, j] = vector[i - 1, j] + mhv[bvindex][j];
                            }
                            else
                            {
                                vector[i, j] = vector[i - 1, j] + mhv[bvindex][j] - N;
                            }
                        }
                    }//for
                }//for
                for (int i = 0; i < N; i++)//N
                {
                    for (int j = 0; j < M - 1; j++)//m-1
                    {
                        vector[i, j] = (2.0 * vector[i, j] - 1) / (2.0 * N);
                    }
                }
                if (bvindex == 0)
                {
                    minindex = bvindex;
                    mincd2 = CalCD2(vector);
                }
                else
                {
                    if (CalCD2(vector) < mincd2)
                    {
                        minindex = bvindex;
                        mincd2 = CalCD2(vector);
                    }
                }
            }//for
            double[,] resultvector;//N*M-1矩阵
            resultvector = new double[N, M - 1];//生成N*M-1 矩阵
            for (int i = 0; i < N; i++)//N
            {
                for (int j = 0; j < M - 1; j++)//M-1
                {
                    if (i == 0)
                    {
                        resultvector[i, j] = mhv[minindex][j];
                    }
                    else
                    {
                        if (resultvector[i - 1, j] + mhv[minindex][j] <= N)
                        {
                            resultvector[i, j] = resultvector[i - 1, j] + mhv[minindex][j];
                        }
                        else
                        {
                            resultvector[i, j] = resultvector[i - 1, j] + mhv[minindex][j] - N ;
                        }
                    }
                }//for
            }//for
            for (int i = 0; i < N; i++)//N
            {
                for (int j = 0; j < M - 1; j++)//m-1
                {
                    resultvector[i, j] = (2.0 * resultvector[i, j] - 1) / (2.0 * N);
                }
            }
            return resultvector;
        }
        public double CalCD2(double[,] m)//计算矩阵的中心化L2-偏差CD2值
        {
            double t1 = Math.Pow((13.0 / 12.0), M - 1);
            double t2 = 0.0;
            double t3 = 0.0;
            for (int k = 0; k < N; k++)
            {
                double temp = 1.0;
                for (int j = 0; j < M - 1; j++)
                {
                    temp *= (1 + 0.5 * (Math.Abs(m[k, j] - 0.5)) - 0.5 * (Math.Pow(Math.Abs(m[k, j] - 0.5), 2)));
                }
                t2 += temp;
            }
            t2 = 2.0 / N * t2;
            for (int k = 0; k < N; k++)
            {
                double t = 0.0;
                for (int j = 0; j < N; j++)
                {
                    double temp = 1.0;
                    for (int i = 0; i < M - 1; i++)
                    {
                        temp *= (1 + 0.5 * (Math.Abs(m[k, i] - 0.5)) + 0.5 * (Math.Abs(m[j, i] - 0.5)) - 0.5 * (Math.Abs(m[k, i] - m[j,i])));
                    }
                    t += temp;
                }
                t3 += t;
            }
            t3 = 1.0 / (N * N)*t3;
            return t1 - t2 + t3; ;
        }
        public double[,] GetMvector()//得到N*M矩阵，这是最后的结果
        {
            double[,] mvector = GetM_1vector();

            double[,] result = new double[N, M];
            for (int i = 0; i < N; i++)
            {
                for (int j = 0; j < M - 1; j++)
                {
                    double cijtemp = 1.0;
                    for (int k = 0; k < j; k++)
                    {
                        cijtemp *= Math.Pow(mvector[i, k], (1.0 / (M - k-1)));
                    }
                    double temp= Math.Pow(mvector[i, j], (1.0 / (M - j-1)));
                    result[i, j] = cijtemp * (1 -temp);
                }
            }
            for (int i = 0; i < N; i++)
            {
                double cijtemp = 1.0;
                for (int k = 0; k < M-1; k++)
                {
                    cijtemp *= Math.Pow(mvector[i, k], 1.0 / (M - k-1));
                }
                result[i, M - 1] = cijtemp;
            }
            return result;
        }
        public int GetNearPrime(int n)//获得与这个数最近的素数
        {
            for (int i = 0, j = 0; ; i++)//, j--)
            {
                if (IsPrime(n + i))
                {
                    return n + i;
                }
                if (IsPrime(n - j))
                {
                    return n - j;
                }
            }
        }
       /// <summary>
        /// 根据权向量的偏好度进行选择权向量
       /// </summary>
       /// <param name="P">子目标的权重,个数当然是M个</param>
       /// <param name="delta">N的变化量N=N*（1+delta）其中0<=delta<=0.4</param>
       /// <returns></returns>
        public double[,] GetMvectorOfP(double[] P,double delta)
        {
            double preference_degree=P[0];
            bool isequal=true;
            for (int i = 0; i < M; i++)
            {
                if (P[i] != preference_degree)
                {
                    isequal=false;
                    break;
                }
            }
            //如果每个子目标的偏好度相等，则直接产生N个
            if (isequal)
            {
                return GetMvector();
            }
            else //如果每个子目标的偏好度不同，则直接产生N*（1+delta）个，然后根据偏好度选取N个
            {
                int PN = N;
                N = (int)(N * (1 + delta));
                double[,] result = new double[PN, M];
                int Nindex=0,Mindex=0;
                double[,] tempresult = new double[N, M];
                tempresult=GetMvector();
                List<int> HavedIndex = new List<int>();
                //挑选权向量的顺序，按照偏好度,Pindex最后里面的结果就是要挑选的子目标的顺序
                int[] Pindex = new int[M];
                for (int i = 0; i < M; i++)
                {
                    Pindex[i] = i;
                }
                    for (int i = 0; i < M; i++)
                    {
                        for (int j = i; j < M; j++)
                        {
                            if (P[i] < P[j])
                            {
                                int tempindex = Pindex[i];
                                Pindex[i] = Pindex[j];
                                Pindex[j] = tempindex;
                            }
                        }
                    }
                //按照偏好度进行选择，选择的顺序就是Pindex里面存储的数据
                    for (int i = 0; i < M; i++)
                    {
                        List<int> MaxIndex = new List<int>();
                        MaxIndex.Clear();
                        MaxIndex = GetMaxNPIndex(tempresult, Pindex[i], (int)(PN * P[Pindex[i]]),ref HavedIndex);
                        Console.WriteLine(MaxIndex.Count+"----------");
                        for (int j = 0; j < MaxIndex.Count; j++)
                        {
                            //将选择的权向量组合进result矩阵中
                            for (int k = 0; k < M; k++)
                            {
                                result[Nindex,Mindex] = tempresult[MaxIndex[j],k];
                                Mindex++;
                            }
                                Mindex = 0;
                                Nindex++;
                                if (Nindex >= N)
                                {
                                    return result;
                                }
                        }                      
                    }
                    return result;
            }
        }
        /// <summary>
        /// 排序，选择每个子目标上的向量index
        /// </summary>
        public List<int> GetMaxNPIndex(double[,] tempresult,int Pindex,int n,ref List<int> havedindex)
        {
            List<int> MaxIndex=new List<int>();
                for (int i = 0; i < N; i++)
                {
                    if (!havedindex.Contains(i))
                    {
                        double max = tempresult[i, Pindex];
                        int index = i;
                        for (int j = i + 1; j < N; j++)
                        {
                            if (!havedindex.Contains(j))
                            {
                                if (tempresult[j, Pindex] > max)
                                {
                                    max = tempresult[j, Pindex];
                                    index = j;
                                }
                            }
                        }
                        MaxIndex.Add(index);
                        havedindex.Add(index);
                        if (MaxIndex.Count >= n)
                        {
                            return MaxIndex;
                        }
                    }
                }
            return null;
        }
    }
}