一种线性约束下的生成随机数修正方法

修正方法

    设需要生成一组如下所示的$N$维随机向量$\mathbf{x}={\left[\begin{array}{ccccc}{x}_1& \cdots & x_j& \cdots & x_N\end{array}\right]}^T$满足如式(1)所示的$N$个随机量上下界约束条件、$M_{eq}$个线性方程约束和$M_{ieq}$个线性不等式约束
{ x j ∈ [ x ‾ j , x ‾ j ] , j ∈ { 1 , ⋯   , N } ∑ j = 1 N a e q , i , j x j = b e q , i , i ∈ { 1 , ⋯   , M e q } ∑ j = 1 N a i e q , i , j x j ≤ b i e q , i , i ∈ { 1 , ⋯   , M i e q } \begin{equation} \begin{cases} x_j ∈ \left[ \underline{x}_j, \overline{x}_j \right], j ∈ \{ 1, \cdots, N \} \\ \displaystyle \sum_{j = 1}^{N} a_{eq, i, j} x_j = b_{eq, i}, i ∈ \{ 1, \cdots, M_{eq} \} \\ \displaystyle \sum_{j = 1}^{N} a_{ieq, i, j} x_j ≤ b_{ieq, i}, i ∈ \{ 1, \cdots, M_{ieq} \} \\ \end{cases} \end{equation} ⎩ ⎨ ⎧​xj​∈[x​j​,xj​],j∈{1,⋯,N}j=1∑N​aeq,i,j​xj​=beq,i​,i∈{1,⋯,Meq​}j=1∑N​aieq,i,j​xj​≤bieq,i​,i∈{1,⋯,Mieq​}​​​
其中：${\underline{x}}_j$和${\overline{x}}_j$分别是第$j$维随机数$x_j$的下界和上界；$a_{eq,i,j}$和$a_{ieq,i,j}$分别为第$i$个线性方程约束和第$i$个线性不等式约束中随机数$x_j$的系数，$a_{eq,i,j}$和$a_{ieq,i,j}$不能全部等于$0$；$b_{eq,i}$和$b_{ieq,i}$分别为第$i$个线性方程和第$i$个线性不等式的常数项。
    若按照如式(2)所示的随机数生成方法进行生成随机数${\mathbf{x}}_{\mathbf{gen}}={\left[\begin{array}{cccc}{x}_{gen,1}& \cdots & x_{gen,i}\cdots & x_{gen,N}\end{array}\right]}^T$作为式(1)所示问题的其中一个测试解
x g e n , j = r j ( x ‾ j − x ‾ j ) + x ‾ j , j ∈ { 1 , ⋯   , N } \begin{equation} x_{gen, j}=r_j \left(\overline{x}_j - \underline{x}_j \right) + \underline{x}_j, j ∈ \{ 1,\cdots, N \} \end{equation} xgen,j​=rj​(xj​−x​j​)+x​j​,j∈{1,⋯,N}​​

其中：$r_j$为$[0,1]$上满足特定概率分布的随机数，易得
x ‾ j ≤ x g e n , j ≤ ( x ‾ j − x ‾ j ) + x ‾ j = x ‾ j , j ∈ { 1 , ⋯   , N } \begin{equation} \underline{x}_j ≤ x_{gen, j} ≤ \left( \overline{x}_j - \underline{x}_j \right) + \underline{x}_j = \underline{x}_j, j ∈ \{ 1, \cdots, N \} \end{equation} x​j​≤xgen,j​≤(xj​−x​j​)+x​j​=x​j​,j∈{1,⋯,N}​​
但该随机数生成方法无法保证满足式(1)中的线性方程约束。当生成的随机向量不满足式(1)中的线性方程约束时，须通过将${\mathbf{x}}_{\mathbf{gen}}$修正为${\mathbf{x}}_{\mathbf{cor}}={\left[\begin{array}{ccccc}{x}_{cor,1}& \cdots & x_{cor,j}& \cdots & x_{cor,N}\end{array}\right]}^T$，使之既满足式(1)中的变量上下界约束，又满足式(1)中的线性方程约束。
    将式(1)中转写成矩阵形式，即有
$$\begin{array}{c}\{\begin{array}{l}\underline{\mathbf{x}}\le \mathbf{x}\le \overline{\mathbf{x}}\\ {\mathbf{A}}_{\mathbf{eq}}\mathbf{x}={\mathbf{b}}_{\mathbf{eq}}\\ {\mathbf{A}}_{\mathbf{ieq}}\mathbf{x}\le {\mathbf{b}}_{\mathbf{ieq}}\end{array}\end{array}$$
其中：矩阵方程中的系数矩阵${\mathbf{A}}_{\mathbf{eq}}=\left[\begin{array}{ccccc}{a}_{eq,1,1}& \cdots & a_{eq,1,j}& \cdots & a_{eq,1,N}\\ ⋮& ⋮& ⋮& ⋮& ⋮\\ a_{eq,i,1}& \cdots & a_{eq,i,j}& \cdots & a_{eq,i,N}\\ ⋮& ⋮& ⋮& ⋮& ⋮\\ a_{eq,M_{eq},1}& \cdots & a_{eq,M_{eq},j}& \cdots & a_{eq,M_{eq},N}\end{array}\right]$，常数向量${\mathbf{b}}_{\mathbf{eq}}=\left[\begin{array}{c}{b}_{eq,1}\\ ⋮\\ b_{eq,i}\\ ⋮\\ b_{eq,M_{eq}}\end{array}\right]$；矩阵不等式中的系数矩阵${\mathbf{A}}_{\mathbf{ieq}}=\left[\begin{array}{ccccc}{a}_{ieq,1,1}& \cdots & a_{ieq,1,j}& \cdots & a_{ieq,1,N}\\ ⋮& ⋮& ⋮& ⋮& ⋮\\ a_{ieq,i,1}& \cdots & a_{ieq,i,j}& \cdots & a_{ieq,i,N}\\ ⋮& ⋮& ⋮& ⋮& ⋮\\ a_{ieq,M_{ieq},1}& \cdots & a_{ieq,M_{ieq},j}& \cdots & a_{ieq,M_{ieq},N}\end{array}\right]$，常数向量${\mathbf{b}}_{\mathbf{ieq}}=\left[\begin{array}{c}{b}_{ieq,1}\\ ⋮\\ b_{ieq,i}\\ ⋮\\ b_{eq,M_{ieq}}\end{array}\right]$。
    将矩阵方程中的系数矩阵${\mathbf{A}}_{\mathbf{eq}}$经过有限次的初等变换即转换成标准形矩阵${\mathbf{A}}_{\mathbf{eq},\mathbf{STD}}$，$\mathbf{x}$通过对应次序的初等列变换（实质作用为调换$\mathbf{x}$中元素的顺序）$\tilde{\mathbf{x}}$，常数向量${\mathbf{b}}_{\mathbf{eq}}$通过对应次序的初等行变换转换成${\mathbf{b}}_{\mathbf{eq},\mathbf{STD}}$。设$R_{eq}\left(R_{eq}\le \min \left(M_{eq},N\right)\right)$为矩阵${\mathbf{A}}_{\mathbf{eq}}$和${\mathbf{A}}_{\mathbf{eq},\mathbf{STD}}$的秩，为简化讨论假定$\mathbf{x}=\tilde{\mathbf{x}}$，则结合线性方程组理论和博客《一种线性方程约束下生成随机数修正的一般性方法（上）》、《一种线性方程约束下生成随机数修正的一般性方法（中）》、《一种线性方程约束下生成随机数修正的一般性方法（下）》中的论述可知，式(1)中线性方程的解有以下情况：

    情况1：若系数向量${\mathbf{b}}_{\mathbf{eq},\mathbf{STD}}$自第$\left(R_{eq}+1\right)$列到第$M_{eq}$列存在非零元素，则式(8)所示的矩阵方程无解，此时无法将${\mathbf{x}}_{\mathbf{gen}}$修正为${\mathbf{x}}_{\mathbf{cor}}$，使之满足式(1)中的线性方程约束。

    情况2：若系数向量${\mathbf{b}}_{\mathbf{eq},\mathbf{STD}}$自第$\left(R_{eq}+1\right)$列到第$M_{eq}$列的元素全为零，则式(8)所示的矩阵方程有解，此时可分为以下情况将${\mathbf{x}}_{\mathbf{gen}}$修正为${\mathbf{x}}_{\mathbf{cor}}$，使得修正后的${\mathbf{x}}_{\mathbf{cor}}$满足式(1)中的线性方程约束：

    情况2.1：若$R_{eq}=1$，此时可将式(4)中约简为
$$\begin{array}{c}\{\begin{array}{l}\underline{\mathbf{x}}\le \mathbf{x}\le \overline{\mathbf{x}}\\ {\mathbf{A}}_{\mathbf{eq},1}\mathbf{x}=b_{eq,1}\\ {\mathbf{A}}_{\mathbf{ieq}}\mathbf{x}\le {\mathbf{b}}_{\mathbf{ieq}}\end{array}\end{array}$$
其中：${\mathbf{A}}_{\mathbf{eq},1}=\left[\begin{array}{ccccc}{a}_{eq,1,1}& \cdots & a_{eq,1,j}& \cdots & a_{eq,1,N}\end{array}\right]$。
    将式(5)进一步转化为如式(6)所示的形式，该式仅包含变量上下界约束和线性不等式约束，然后利用博客《一种线性不等式约束下的生成随机数修正方法》所述方法即可将${\mathbf{x}}_{\mathbf{gen}}$修正为${\mathbf{x}}_{\mathbf{cor}}$。
$$\begin{array}{c}\{\begin{array}{l}\underline{\mathbf{x}}\le \mathbf{x}\le \overline{\mathbf{x}}\\ {\mathbf{A}}_{\mathbf{eq},1}\mathbf{x}\le b_{eq,1}\\ -{\mathbf{A}}_{\mathbf{eq},1}\mathbf{x}\le -b_{eq,1}\\ {\mathbf{A}}_{\mathbf{ieq}}\mathbf{x}\le {\mathbf{b}}_{\mathbf{ieq}}\end{array}\end{array}$$

    情况2.2：若$R_{eq}>1$，则对于式(4)中的矩阵方程部分，根据博客《一种线性方程约束下生成随机数修正的一般性方法（上）》中的论述可知，${\mathbf{A}}_{\mathbf{eq},\mathbf{STD}}$的结构如式(7)所示
$$\begin{array}{c}{\mathbf{A}}_{\mathbf{eq},\mathbf{STD}}=\left[\begin{array}{cc}{\mathbf{E}}_{R_{eq}}& {{\tilde{\mathbf{A}}}_{\mathbf{eq}}}_{R_{eq}\times \left(N-R_{eq}\right)}\\ {\mathbf{O}}_{\left(M_{eq}-R_{eq}\right)\times R_{eq}}& {\mathbf{O}}_{\left(M_{eq}-R_{eq}\right)\times \left(N-R_{eq}\right)}\end{array}\right]\end{array}$$
其中：${\mathbf{E}}_{R_{eq}}$为$R_{eq}$阶单位矩阵，${{\tilde{\mathbf{A}}}_{\mathbf{eq}}}_{R_{eq}\times \left(N-R_{eq}\right)}=\left[\begin{array}{ccccc}{\tilde{a}}_{eq,1,1}& \cdots & {\tilde{a}}_{eq,1,j}& \cdots & {\tilde{a}}_{eq,1,\left(N-R_{eq}\right)}\\ ⋮& ⋮& ⋮& ⋮& ⋮\\ {\tilde{a}}_{eq,i,1}& \cdots & {\tilde{a}}_{eq,i,j}& \cdots & {\tilde{a}}_{eq,i,\left(N-R_{eq}\right)}\\ ⋮& ⋮& ⋮& ⋮& ⋮\\ {\tilde{a}}_{eq,R_{eq},1}& \cdots & {\tilde{a}}_{eq,R_{eq},j}& \cdots & {\tilde{a}}_{eq,R_{eq},\left(N-R_{eq}\right)}\end{array}\right]$，${\mathbf{O}}_{\left(M_{eq}-R_{eq}\right)\times R_{eq}}$和${\mathbf{O}}_{\left(M_{eq}-R_{eq}\right)\times \left(N-R_{eq}\right)}$分别为$\left(M_{eq}-R_{eq}\right)$行$R_{eq}$列和$\left(M_{eq}-R_{eq}\right)$行$\left(N-R_{eq}\right)$列的零矩阵。
    进一步地，可将式(4)中的矩阵方程部分转换为如式(7)所示的矩阵不等式
$$\begin{array}{c}\{\begin{array}{l}\underline{\mathbf{x}}\le \mathbf{x}\le \overline{\mathbf{x}}\\ {\mathbf{b}}_{\mathbf{eq},\mathbf{STD}}-{\overline{\mathbf{x}}}_{1-{\mathbf{R}}_{\mathbf{eq}}}\le {\tilde{\mathbf{A}}}_{\mathbf{eq}}{\mathbf{x}}_{1-{\mathbf{R}}_{\mathbf{eq}}}\le {\mathbf{b}}_{\mathbf{eq},\mathbf{STD}}-{\underline{\mathbf{x}}}_{1-{\mathbf{R}}_{\mathbf{eq}}}\end{array}\end{array}$$
其中：${\mathbf{x}}_{1-{\mathbf{R}}_{\mathbf{eq}}}={\left[\begin{array}{ccccc}{x}_1& \cdots & x_j& \cdots & x_{R_{eq}}\end{array}\right]}^T$，${\underline{\mathbf{x}}}_{1-{\mathbf{R}}_{\mathbf{eq}}}={\left[\begin{array}{ccccc}{\underline{x}}_1& \cdots & {\underline{x}}_j& \cdots & {\underline{x}}_{R_{eq}}\end{array}\right]}^T$，${\overline{\mathbf{x}}}_{1-{\mathbf{R}}_{\mathbf{eq}}}={\left[\begin{array}{ccccc}{\overline{x}}_1& \cdots & {\overline{x}}_j& \cdots & {\overline{x}}_{R_{eq}}\end{array}\right]}^T$。
    对于式(4)中的矩阵不等式部分，将其进行矩阵的分块，并在此基础上进一步转化可得如式(9)所示形式
$$\begin{array}{c}\begin{array}{rl}{\mathbf{A}}_{\mathbf{ieq}}\mathbf{x}& =\left[\begin{array}{cc}{\mathbf{A}}_{\mathbf{ieq},1-{\mathbf{R}}_{\mathbf{eq}}}& {\mathbf{A}}_{\mathbf{ieq},\left({\mathbf{R}}_{\mathbf{eq}}+1\right)-\mathbf{N}}\end{array}\right]\left[\begin{array}{c}{\mathbf{x}}_{1-{\mathbf{R}}_{\mathbf{eq}}}\\ {\mathbf{x}}_{\left({\mathbf{R}}_{\mathbf{eq}}+1\right)-\mathbf{N}}\end{array}\right]\\ & ={\mathbf{A}}_{\mathbf{ieq},1-{\mathbf{R}}_{\mathbf{eq}}}{\mathbf{x}}_{1-{\mathbf{R}}_{\mathbf{eq}}}+{\mathbf{A}}_{\mathbf{ieq},\left({\mathbf{R}}_{\mathbf{eq}}+1\right)-\mathbf{N}}{\mathbf{x}}_{\left({\mathbf{R}}_{\mathbf{eq}}+1\right)-\mathbf{N}}\\ & ={\mathbf{A}}_{\mathbf{ieq},1-{\mathbf{R}}_{\mathbf{eq}}}\left({\mathbf{b}}_{\mathbf{eq},\mathbf{STD}}-{\tilde{\mathbf{A}}}_{\mathbf{eq}}{\mathbf{x}}_{\left({\mathbf{R}}_{\mathbf{eq}}+1\right)-\mathbf{N}}\right)+{\mathbf{A}}_{\mathbf{ieq},\left({\mathbf{R}}_{\mathbf{eq}}+1\right)-\mathbf{N}}{\mathbf{x}}_{\left({\mathbf{R}}_{\mathbf{eq}}+1\right)-\mathbf{N}}\\ & ={\mathbf{A}}_{\mathbf{ieq},1-{\mathbf{R}}_{\mathbf{eq}}}{\mathbf{b}}_{\mathbf{eq},\mathbf{STD}}+\left({\mathbf{A}}_{\mathbf{ieq},\left({\mathbf{R}}_{\mathbf{eq}}+1\right)-\mathbf{N}}-{\mathbf{A}}_{\mathbf{ieq},1-{\mathbf{R}}_{\mathbf{eq}}}{\tilde{\mathbf{A}}}_{\mathbf{eq}}\right){\mathbf{x}}_{\left({\mathbf{R}}_{\mathbf{eq}}+1\right)-\mathbf{N}}\\ & \le {\mathbf{b}}_{\mathbf{ieq}}\end{array}\end{array}$$
其中：${\mathbf{A}}_{\mathbf{ieq},1-{\mathbf{R}}_{\mathbf{eq}}}=\left[\begin{array}{ccccc}{a}_{ieq,1,1}& \cdots & a_{ieq,1,j}& \cdots & a_{ieq,1,R_{eq}}\\ ⋮& ⋮& ⋮& ⋮& ⋮\\ a_{ieq,i,1}& \cdots & a_{ieq,i,j}& \cdots & a_{ieq,i,R_{eq}}\\ ⋮& ⋮& ⋮& ⋮& ⋮\\ a_{ieq,M_{ieq},1}& \cdots & a_{ieq,M_{ieq},j}& \cdots & a_{ieq,M_{ieq},R_{eq}}\end{array}\right]$，${\mathbf{A}}_{\mathbf{ieq},\left({\mathbf{R}}_{\mathbf{eq}}+1\right)-\mathbf{N}}=\left[\begin{array}{ccccc}{a}_{ieq,1,R_{eq}+1}& \cdots & a_{ieq,1,j}& \cdots & a_{ieq,1,N}\\ ⋮& ⋮& ⋮& ⋮& ⋮\\ a_{ieq,i,R_{eq}+1}& \cdots & a_{ieq,i,j}& \cdots & a_{ieq,i,N}\\ ⋮& ⋮& ⋮& ⋮& ⋮\\ a_{ieq,M_{ieq},R_{eq}+1}& \cdots & a_{ieq,M_{ieq},j}& \cdots & a_{ieq,M_{ieq},N}\end{array}\right]$；${\mathbf{x}}_{\left({\mathbf{R}}_{\mathbf{eq}}+1\right)-\mathbf{N}}={\left[\begin{array}{ccccc}{x}_{R_{eq}+1}& \cdots & x_j& \cdots & x_N\end{array}\right]}^T$，${\underline{\mathbf{x}}}_{\left({\mathbf{R}}_{\mathbf{eq}}+1\right)-\mathbf{N}}={\left[\begin{array}{ccccc}{\underline{x}}_{R_{eq}+1}& \cdots & {\underline{x}}_j& \cdots & {\underline{x}}_N\end{array}\right]}^T$，${\overline{\mathbf{x}}}_{\left({\mathbf{R}}_{\mathbf{eq}}+1\right)-\mathbf{N}}={\left[\begin{array}{ccccc}{\overline{x}}_{R_{eq}+1}& \cdots & {\overline{x}}_j& \cdots & {\overline{x}}_N\end{array}\right]}^T$，${\mathbf{x}}_{1-{\mathbf{R}}_{\mathbf{eq}}}$与${\underline{\mathbf{x}}}_{\left({\mathbf{R}}_{\mathbf{eq}}+1\right)-\mathbf{N}}$之间满足
$$\begin{array}{c}{\mathbf{x}}_{1-{\mathbf{R}}_{\mathbf{eq}}}={\mathbf{b}}_{\mathbf{eq},\mathbf{STD}}-{\tilde{\mathbf{A}}}_{\mathbf{eq}}{\mathbf{x}}_{\left({\mathbf{R}}_{\mathbf{eq}}+1\right)-\mathbf{N}}\end{array}$$
    对式(9)进行移项可得
$$\begin{array}{c}\left({\mathbf{A}}_{\mathbf{ieq},\left({\mathbf{R}}_{\mathbf{eq}}+1\right)-\mathbf{N}}-{\mathbf{A}}_{\mathbf{ieq},1-{\mathbf{R}}_{\mathbf{eq}}}{\tilde{\mathbf{A}}}_{\mathbf{eq}}\right){\mathbf{x}}_{\left({\mathbf{R}}_{\mathbf{eq}}+1\right)-\mathbf{N}}\le {\mathbf{b}}_{\mathbf{ieq}}-{\mathbf{A}}_{\mathbf{ieq},1-{\mathbf{R}}_{\mathbf{eq}}}{\mathbf{b}}_{\mathbf{eq},\mathbf{STD}}\end{array}$$
将其与式(8)联立后，再进行如式(12)所示的进一步转化即可变为以${\mathbf{x}}_{\left({\mathbf{R}}_{\mathbf{eq}}+1\right)-\mathbf{N}}$为变量的矩阵不等式
$$\begin{array}{c}\{\begin{array}{l}{\underline{\mathbf{x}}}_{\left({\mathbf{R}}_{\mathbf{eq}}+1\right)-\mathbf{N}}\le {\mathbf{x}}_{\left({\mathbf{R}}_{\mathbf{eq}}+1\right)-\mathbf{N}}\le {\overline{\mathbf{x}}}_{\left({\mathbf{R}}_{\mathbf{eq}}+1\right)-\mathbf{N}}\\ -{\tilde{\mathbf{A}}}_{\mathbf{eq}}{\mathbf{x}}_{\left({\mathbf{R}}_{\mathbf{eq}}+1\right)-\mathbf{N}}\le -\left({\mathbf{b}}_{\mathbf{eq},\mathbf{STD}}-{\overline{\mathbf{x}}}_{1-{\mathbf{R}}_{\mathbf{eq}}}\right)\\ {\tilde{\mathbf{A}}}_{\mathbf{eq}}{\mathbf{x}}_{\left({\mathbf{R}}_{\mathbf{eq}}+1\right)-\mathbf{N}}\le {\mathbf{b}}_{\mathbf{eq},\mathbf{STD}}-{\underline{\mathbf{x}}}_{1-{\mathbf{R}}_{\mathbf{eq}}}\\ \left({\mathbf{A}}_{\mathbf{ieq},\left({\mathbf{R}}_{\mathbf{eq}}+1\right)-\mathbf{N}}-{\mathbf{A}}_{\mathbf{ieq},1-{\mathbf{R}}_{\mathbf{eq}}}{\tilde{\mathbf{A}}}_{\mathbf{eq}}\right){\mathbf{x}}_{\left({\mathbf{R}}_{\mathbf{eq}}+1\right)-\mathbf{N}}\le {\mathbf{b}}_{\mathbf{ieq}}-{\mathbf{A}}_{\mathbf{ieq},1-{\mathbf{R}}_{\mathbf{eq}}}{\mathbf{b}}_{\mathbf{eq},\mathbf{STD}}\end{array}\end{array}$$
随后利用博客《一种线性不等式约束下的生成随机数修正方法》所述方法将$x_{gen,R_{eq}+1},\cdots ,x_{gen,N}$修正为$x_{cor,R_{eq}+1},\cdots ,x_{cor,N}$，再根据式(10)将$x_{gen,1},\cdots ,x_{gen,R_{eq}}$修正为$x_{cor,1},\cdots ,x_{cor,R_{eq}}$，从而将${\mathbf{x}}_{\mathbf{gen}}$修正为${\mathbf{x}}_{\mathbf{cor}}$。

MATLAB代码

linear_inequalities_constraints_correction函数

function [X_cor, correction_info] = linear_inequalities_constraints_correction(X_gen, lb, ub, A_ieq, b_ieq, ...
                                    W_multiplicator)
    % LINEAR_INEQUALITIES_CONSTRAINTS_CORRECTION  Correction onto generated random vector.
    %                                             A methodology of correction onto generated random vector,
    %                                             facilitating the corrected vector to range within variable
    %                                             boundaries and to meet linear constraints

    if ~exist('W_multiplicator', 'var')
        W_multiplicator = 1000;
    end
    max_runtime = ceil(W_multiplicator / size(X_gen, 2));

    M = size(A_ieq, 1);
    [N_D, N_P] = size(X_gen);
    
    b_ieq_lb = diag(A_ieq * ((A_ieq > 0)' .* lb + (A_ieq < 0)' .* ub));
    b_ieq_ub = min([diag(A_ieq * ((A_ieq > 0)' .* ub + (A_ieq < 0)' .* lb)), b_ieq], [], 2);
    
    if any(b_ieq < b_ieq_lb)
        X_cor = [];
        correction_info.is_feasible = false(1, N_P);
        correction_info.does_X_cor_exceed.by_bit = true(N_D, N_P);
        correction_info.does_X_cor_exceed.globally = true(1, N_P);
        correction_info.is_inequality_less_than_b.by_bit = false(M, N_P);
        correction_info.is_inequality_less_than_b.globally = false;

    else
        X_cor = X_gen;

        run = 0;
        while 1
            X_cor = X_gen;
            iEq_gen = A_ieq * X_gen;
            does_A_ieq_X_gen_exceed = (iEq_gen > b_ieq);
    
            for indi = 1:N_P
                if any(does_A_ieq_X_gen_exceed(:, indi))
                    flag_A_ieq_X_gen_exceed = find(does_A_ieq_X_gen_exceed(:, indi));
                    row_A_ieq_candi = flag_A_ieq_X_gen_exceed(randi(length(flag_A_ieq_X_gen_exceed)));
                    A_ieq_candi = A_ieq(row_A_ieq_candi, :);
                    Eq_candi_gen = A_ieq_candi * X_gen(:, indi);
                    b_ieq_k_lb = b_ieq_lb(row_A_ieq_candi, :);
                    b_ieq_k_ub = b_ieq_ub(row_A_ieq_candi, :);
                    delta_plus = A_ieq_candi * ((A_ieq_candi > 0)' .* (X_gen(:, indi) - lb) ...
                        + (A_ieq_candi < 0)' .* (X_gen(:, indi) - ub));
                    x_cor_base = zeros(N_D, 1);
    
                    if Eq_candi_gen > b_ieq_k_ub
                        b_eq_k_best = [b_ieq_k_lb + (delta_plus * (b_ieq_lb - A_ieq * ((A_ieq_candi > 0)' ...
                            .* lb + (A_ieq_candi < 0)' .* ub + (A_ieq_candi == 0)' .* X_gen(:, indi)))) ...
                            ./ (A_ieq * ((A_ieq_candi > 0)' .* (X_gen(:, indi) - lb) ...
                            + (A_ieq_candi < 0)' .* (X_gen(:, indi) - ub))), ...
                            b_ieq_k_lb + (delta_plus * (b_ieq_ub - A_ieq * ((A_ieq_candi > 0)' ...
                            .* lb + (A_ieq_candi < 0)' .* ub + (A_ieq_candi == 0)' .* X_gen(:, indi)))) ...
                            ./ (A_ieq * ((A_ieq_candi > 0)' .* (X_gen(:, indi) - lb) ...
                            + (A_ieq_candi < 0)' .* (X_gen(:, indi) - ub)))];
                        b_eq_k_min = max(min(b_eq_k_best, [], 2), [], 1);
                        b_eq_k_max = min(max(b_eq_k_best, [], 2), [], 1);
                        if b_eq_k_min <= b_eq_k_max
                            b_eq_k = b_eq_k_min + rand * (b_eq_k_max - b_eq_k_min);
                        else
                            b_eq_k = b_eq_k_best(row_A_ieq_candi, 1) + rand * (b_eq_k_best(row_A_ieq_candi, 2) ...
                                     - b_eq_k_best(row_A_ieq_candi, 1));
                        end
                        x_cor_base(A_ieq_candi > 0) = (((b_eq_k - b_ieq_k_lb) / delta_plus) ...
                                                      * (X_gen(A_ieq_candi > 0, indi) - lb(A_ieq_candi > 0)) ...
                                                      + lb(A_ieq_candi > 0))';
                        x_cor_base(A_ieq_candi == 0) = (X_gen(A_ieq_candi == 0))';
                        x_cor_base(A_ieq_candi < 0) = (((b_eq_k - b_ieq_k_lb) / delta_plus) ...
                                                      * (X_gen(A_ieq_candi < 0, indi) - ub(A_ieq_candi < 0)) ...
                                                      + ub(A_ieq_candi < 0))';
                    end
                    X_cor(:, indi) = x_cor_base;
                end
            end
    
            for indi = 1:N_P
                X_cor(abs(X_cor(:, indi) - lb) <= 1e-12, indi) = lb(abs(X_cor(:, indi) - lb) <= 1e-12);
                X_cor(abs(X_cor(:, indi) - ub) <= 1e-12, indi) = ub(abs(X_cor(:, indi) - ub) <= 1e-12);
            end
    
            correction_info.is_feasible = all((X_cor - lb >= -1e-12) & (X_cor - ub <= 1e-12)) ...
                                          & all((A_ieq * X_cor - b_ieq) <= 1e-12);
            correction_info.does_X_cor_exceed.by_bit = (X_cor - lb < -1e-12) | (X_cor - ub > 1e-12);
            correction_info.does_X_cor_exceed.globally = any((X_cor - lb < -1e-12) | (X_cor - ub > 1e-12));
            correction_info.is_inequality_less_than_b.by_bit = (A_ieq * X_cor - b_ieq) <= 1e-12;
            correction_info.is_inequality_less_than_b.globally = all((A_ieq * X_cor - b_ieq) <= 1e-12);

            run = run + 1;
            % Feasible, Theorem 5
            if all(correction_info.is_feasible)
                break;
            end
            % Infeasible, Theorems 4 & 5
            if ~any(correction_info.is_feasible) && run >= max_runtime
                X_cor = [];
                correction_info.is_feasible = false(1, N_P);
                correction_info.does_X_cor_exceed.by_bit = true(N_D, N_P);
                correction_info.does_X_cor_exceed.globally = true(1, N_P);
                correction_info.is_equation_equal_to_b.by_bit = false(M, N_P);
                correction_info.is_equation_equal_to_b.globally = false;
                break;
            end

            X_gen = X_cor; 
        end

    end

end

linear_constraints_correction函数

function [X_cor, correction_info] = linear_constraints_correction(X_gen, lb, ub, A_eq, b_eq, A_ieq, b_ieq, ...
                                    W_multiplicator)
    % LINEAR_PROGRAMMING_CORRECTION  Correction onto generated random vector.
    %                                A methodology of correction onto generated random vector,
    %                                facilitating the corrected vector to range within variable
    %                                boundaries and to meet linear constraints

    if ~exist('W_multiplicator', 'var')
        W_multiplicator = 1000;
    end

    Ab_eq = [A_eq, b_eq];
    
    M_eq = size(A_eq, 1);
    M_ieq = size(A_ieq, 1);
    [N_D, N_P] = size(X_gen);

    % Invoke Python for an accurate computation of the reduced row echelon form of Ab_eq
    Ab_eq_py = py.numpy.array(Ab_eq);
    Ab_eq_RREF_sym = py.sympy.Matrix(Ab_eq_py).rref();
    Ab_eq_RREF_raw = double(py.numpy.asarray(Ab_eq_RREF_sym(1), dtype = 'float'));
    if size(Ab_eq, 1) == 1
        Ab_eq_RREF = Ab_eq_RREF_raw;
    else
        Ab_eq_RREF = reshape(Ab_eq_RREF_raw, size(Ab_eq_RREF_raw, 2), size(Ab_eq_RREF_raw, 3));
    end

    R_A_eq = sum(any(Ab_eq_RREF(:, 1:end - 1), 2), 1);
    R_Ab_eq = sum(any(Ab_eq_RREF, 2), 1);

    % pivots_eq: Nonzero pivot columns in Ab_eq_RREF
    % non_pivots_eq: Other columns, except the end column, in Ab_eq_RREF
    pivots_eq = zeros(1, R_A_eq);
    for count = 1:R_A_eq
        pivots_eq(1, count) = find(Ab_eq_RREF(count, :), 1);
    end
    non_pivots_eq = find(~ismember(1:N_D, pivots_eq));
    
    A_eq_tilde = Ab_eq_RREF(:, non_pivots_eq);
    b_eq_tilde = Ab_eq_RREF(:, end);
    
    if R_A_eq < R_Ab_eq
        X_cor = [];
        correction_info.is_feasible = false(1, N_P);
        correction_info.does_X_cor_exceed.by_bit = true(N_D, N_P);
        correction_info.does_X_cor_exceed.globally = true(1, N_P);
        correction_info.is_equation_equal_to_b.by_bit = false(M_eq, N_P);
        correction_info.is_equation_equal_to_b.globally = false;
        correction_info.is_inequality_less_than_b.by_bit = false(M_ieq, N_P);
        correction_info.is_inequality_less_than_b.globally = false;
    
    else
        lb_base = lb(non_pivots_eq, :);
        ub_base = ub(non_pivots_eq, :);
        b_eq_tilde_plus = b_eq_tilde - lb(pivots_eq, end);
        b_eq_tilde_minus = b_eq_tilde - ub(pivots_eq, end);
        if isempty(A_eq)
            A_new = A_ieq;
            b_new = b_ieq;
            [X_cor, ~] = linear_inequalities_constraints_correction(X_gen, lb, ub, A_new, ...
                         b_new, W_multiplicator);

        elseif isempty(A_ieq)
            if R_A_eq == 1
                A_new = [-A_eq; A_eq];
                b_new = [-b_eq; b_eq];
                [X_cor, ~] = linear_inequalities_constraints_correction(X_gen, lb, ub, A_new, ...
                             b_new, W_multiplicator);

            else
                A_new = [-A_eq_tilde; A_eq_tilde];
                b_new = [-b_eq_tilde_minus; b_eq_tilde_plus];
                X_gen_base = X_gen(non_pivots_eq, :);
                [X_cor_base, ~] = linear_inequalities_constraints_correction(X_gen_base, lb_base, ...
                                  ub_base, A_new, b_new, W_multiplicator);
                X_cor(pivots_eq, :) = b_eq_tilde - A_eq_tilde * X_cor_base;
                X_cor(non_pivots_eq, :) = X_cor_base;

            end

        else
            if R_A_eq == 1
                A_new = [-A_eq; A_eq; A_ieq];
                b_new = [-b_eq; b_eq; b_ieq];
                [X_cor, ~] = linear_inequalities_constraints_correction(X_gen, lb, ub, A_new, ...
                             b_new, W_multiplicator);

            else
                A_new = [-A_eq_tilde; A_eq_tilde; A_ieq(:, non_pivots_eq) - A_ieq(:, pivots_eq) * A_eq_tilde];
                b_new = [-b_eq_tilde_minus; b_eq_tilde_plus; b_ieq - A_ieq(:, pivots_eq) * b_eq_tilde];
                X_gen_base = X_gen(non_pivots_eq, :);
                [X_cor_base, ~] = linear_inequalities_constraints_correction(X_gen_base, lb_base, ...
                                  ub_base, A_new, b_new, W_multiplicator);
                X_cor(pivots_eq, :) = b_eq_tilde - A_eq_tilde * X_cor_base;
                X_cor(non_pivots_eq, :) = X_cor_base;

            end
        end
        
        if isempty(A_eq)
            correction_info.is_feasible = all((X_cor - lb >= -1e-12) & (X_cor - ub <= 1e-12)) ...
                                          & all((A_ieq * X_cor - b_ieq) <= 1e-12);
            correction_info.does_X_cor_exceed.by_bit = (X_cor - lb < -1e-12) | (X_cor - ub > 1e-12);
            correction_info.does_X_cor_exceed.globally = any((X_cor - lb < -1e-12) | (X_cor - ub > 1e-12));
            correction_info.is_equation_equal_to_b.by_bit = [];
            correction_info.is_equation_equal_to_b.globally = [];
            correction_info.is_inequality_less_than_b.by_bit = (A_ieq * X_cor - b_ieq) <= 1e-12;
            correction_info.is_inequality_less_than_b.globally = all((A_ieq * X_cor - b_ieq) <= 1e-12);

        elseif isempty(A_ieq)
            correction_info.is_feasible = all((X_cor - lb >= -1e-12) & (X_cor - ub <= 1e-12)) ...
                                          & all(abs(A_eq * X_cor - b_eq) <= 1e-12);
            correction_info.does_X_cor_exceed.by_bit = (X_cor - lb < -1e-12) | (X_cor - ub > 1e-12);
            correction_info.does_X_cor_exceed.globally = any((X_cor - lb < -1e-12) | (X_cor - ub > 1e-12));
            correction_info.is_equation_equal_to_b.by_bit = abs(A_eq * X_cor - b_eq) <= 1e-12;
            correction_info.is_equation_equal_to_b.globally = all(abs(A_eq * X_cor - b_eq) <= 1e-12);
            correction_info.is_inequality_less_than_b.by_bit = [];
            correction_info.is_inequality_less_than_b.globally = [];

        else
            correction_info.is_feasible = all((X_cor - lb >= -1e-12) & (X_cor - ub <= 1e-12)) ...
                                          & all(abs(A_eq * X_cor - b_eq) <= 1e-12) & ...
                                          all((A_ieq * X_cor - b_ieq) <= 1e-12);
            correction_info.does_X_cor_exceed.by_bit = (X_cor - lb < -1e-12) | (X_cor - ub > 1e-12);
            correction_info.does_X_cor_exceed.globally = any((X_cor - lb < -1e-12) | (X_cor - ub > 1e-12));
            correction_info.is_equation_equal_to_b.by_bit = abs(A_eq * X_cor - b_eq) <= 1e-12;
            correction_info.is_equation_equal_to_b.globally = all(abs(A_eq * X_cor - b_eq) <= 1e-12);
            correction_info.is_inequality_less_than_b.by_bit = (A_ieq * X_cor - b_ieq) <= 1e-12;
            correction_info.is_inequality_less_than_b.globally = all((A_ieq * X_cor - b_ieq) <= 1e-12);

        end

    end

end

Problem 1

clear;
close all;
clc;

lb = [-3, -1, 2, -4, 0, 1, -6, -5]';
ub = [6, 4, 10, 8, 11, 3, 4, 4]';
X_gen = lb + rand(1, 100) .* (ub - lb);
A_eq = [6, -3, 1, 7, -5, 0, -3, 8];
b_eq = 12;

[X_cor, correction_info] = linear_constraints_correction(X_gen, lb, ub, [], [], A_eq, b_eq);

Problem 2

clear;
close all;
clc;

lb = [-5, -3, -6, 0]';
ub = [7, 1, -2, 6]';
X_gen = lb + rand(1, 100) .* (ub - lb);
A_eq = [1, 1, 1, 1; 0, 2, 1, 1];
b_eq = [-1, 1]';

[X_cor, correction_info] = linear_constraints_correction(X_gen, lb, ub, A_eq, b_eq, [], []);

Problem 3

clear;
close all;
clc;

lb = [-8, -15, -2, 0, -3, -10]';
ub = [9, 7, 4, 5, 8, 2]';
X_gen = lb + rand(1, 100) .* (ub - lb);
A_eq = [4, 3, 5, -7, 6, -8];
b_eq = -2;
A_ieq = [-7, -4, 8, -5, 0, 1];
b_ieq = 14;

[X_cor, correction_info] = linear_constraints_correction(X_gen, lb, ub, A_eq, b_eq, A_ieq, b_ieq);

Problem 4

clear;
close all;
clc;

lb = [-8, -15, -2, 1, -3, -10]';
ub = [9, 7, 11, 7, 8, 2]';
X_gen = lb + rand(1, 100) .* (ub - lb);
A_eq = [4, 3, 5, -7, 6, -8; -7, -4, 8, -5, 0, 8];
b_eq = [-2; 14];
A_ieq = [10, 3, -3, 6, -7, 2; 2, -3, -7, 5, -6, 3];
b_ieq = [9; -8];

[X_cor, correction_info] = linear_constraints_correction(X_gen, lb, ub, A_eq, b_eq, A_ieq, b_ieq);

研究目标

    (1) 探究满足某些非线性约束的生成随机数修正方法；
    (2) 探究随机向量为离散向量和连续-离散混合向量时的生成随机数修正方法；
    (3) 将生成随机数修正方法应用到基于启发式算法的优化问题求解中，使得启发式算法能够始终在优化问题的可行域中搜寻问题的解，从而加快启发式优化算法的收敛速度。