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

    本文内容接续博客 《一种线性方程约束下生成随机数修正的一般性方法（中）》

修正方法

    情况2.2.2.2.2.2.1：当$\sum _{j=1,{\tilde{a}}_{k,j}>0}^{N-R}{\tilde{a}}_{k,j}\left(x_{gen,R+j}-{\overline{x}}_{R+j}\right)+\sum _{j=1,{\tilde{a}}_{k,j}<0}^{N-R}{\tilde{a}}_{k,j}\left(x_{gen,R+j}-{\underline{x}}_{R+j}\right)$与$\sum _{j=1,{\tilde{a}}_{i,j}>0}^{N-R}{\tilde{a}}_{i,j}\left(x_{gen,R+j}-{\overline{x}}_{R+j}\right)+\sum _{j=1,{\tilde{a}}_{i,j}<0}^{N-R}{\tilde{a}}_{i,j}\left(x_{gen,R+j}-{\underline{x}}_{R+j}\right)$同号时，${\underline{\tilde{b}}}_{cor,k}^{\prime }\left(i\right)$、${\overline{\tilde{b}}}_{cor,k}^{\prime }\left(i\right)$、${\underline{\tilde{b}}}_{cor,k}^{\prime \prime }\left(i\right)$、${\overline{\tilde{b}}}_{cor,k}^{\prime \prime }\left(i\right)$分别满足
$$\begin{array}{c}\begin{array}{rl}{\underline{\tilde{b}}}_{cor,k}^{\prime }\left(i\right)& =\left(\sum _{j=1,{\tilde{a}}_{k,j}>0}^{N-R}{\tilde{a}}_{i,j}{\overline{x}}_{R+j}+\sum _{j=1,{\tilde{a}}_{k,j}<0}^{N-R}{\tilde{a}}_{i,j}{\underline{x}}_{R+j}\right)+\frac{\sum _{j=1,{\tilde{a}}_{k,j}>0}^{N-R}{\tilde{a}}_{k,j}\left(x_{gen,R+j}-{\underline{x}}_{R+j}\right)+\sum _{j=1,{\tilde{a}}_{k,j}<0}^{N-R}{\tilde{a}}_{k,j}\left(x_{gen,R+j}-{\overline{x}}_{R+j}\right)}{\sum _{j=1,{\tilde{a}}_{k,j}>0}^{N-R}{\tilde{a}}_{i,j}\left(x_{gen,R+j}-{\underline{x}}_{R+j}\right)+\sum _{j=1,{\tilde{a}}_{k,j}<0}^{N-R}{\tilde{a}}_{i,j}\left(x_{gen,R+j}-{\overline{x}}_{R+j}\right)}\times \left[\left(\sum _{j=1,{\tilde{a}}_{i,j}>0}^{N-R}{\tilde{a}}_{i,j}{\underline{x}}_{R+j}+\sum _{j=1,{\tilde{a}}_{i,j}<0}^{N-R}{\tilde{a}}_{i,j}{\overline{x}}_{R+j}\right)-\left(\sum _{j=1,{\tilde{a}}_{k,j}>0}^{N-R}{\tilde{a}}_{i,j}{\underline{x}}_{R+j}+\sum _{j=1,{\tilde{a}}_{k,j}<0}^{N-R}{\tilde{a}}_{i,j}{\overline{x}}_{R+j}+\sum _{j=1,{\tilde{a}}_{k,j}=0}^{N-R}{\tilde{a}}_{i,j}{x}_{gen,R+j}\right)\right]\\ {\overline{\tilde{b}}}_{cor,k}^{\prime }\left(i\right)& =\left(\sum _{j=1,{\tilde{a}}_{k,j}>0}^{N-R}{\tilde{a}}_{i,j}{\overline{x}}_{R+j}+\sum _{j=1,{\tilde{a}}_{k,j}<0}^{N-R}{\tilde{a}}_{i,j}{\underline{x}}_{R+j}\right)+\frac{\sum _{j=1,{\tilde{a}}_{k,j}>0}^{N-R}{\tilde{a}}_{k,j}\left(x_{gen,R+j}-{\underline{x}}_{R+j}\right)+\sum _{j=1,{\tilde{a}}_{k,j}<0}^{N-R}{\tilde{a}}_{k,j}\left(x_{gen,R+j}-{\overline{x}}_{R+j}\right)}{\sum _{j=1,{\tilde{a}}_{k,j}>0}^{N-R}{\tilde{a}}_{i,j}\left(x_{gen,R+j}-{\underline{x}}_{R+j}\right)+\sum _{j=1,{\tilde{a}}_{k,j}<0}^{N-R}{\tilde{a}}_{i,j}\left(x_{gen,R+j}-{\overline{x}}_{R+j}\right)}\times \left[\left(\sum _{j=1,{\tilde{a}}_{i,j}>0}^{N-R}{\tilde{a}}_{i,j}{\overline{x}}_{R+j}+\sum _{j=1,{\tilde{a}}_{i,j}<0}^{N-R}{\tilde{a}}_{i,j}{\underline{x}}_{R+j}\right)-\left(\sum _{j=1,{\tilde{a}}_{k,j}>0}^{N-R}{\tilde{a}}_{i,j}{\underline{x}}_{R+j}+\sum _{j=1,{\tilde{a}}_{k,j}<0}^{N-R}{\tilde{a}}_{i,j}{\overline{x}}_{R+j}+\sum _{j=1,{\tilde{a}}_{k,j}=0}^{N-R}{\tilde{a}}_{i,j}{x}_{gen,R+j}\right)\right]\\ {\underline{\tilde{b}}}_{cor,k}^{\prime \prime }\left(i\right)& =\left(\sum _{j=1,{\tilde{a}}_{k,j}>0}^{N-R}{\tilde{a}}_{i,j}{\overline{x}}_{R+j}+\sum _{j=1,{\tilde{a}}_{k,j}<0}^{N-R}{\tilde{a}}_{i,j}{\underline{x}}_{R+j}\right)+\frac{\sum _{j=1,{\tilde{a}}_{k,j}>0}^{N-R}{\tilde{a}}_{k,j}\left(x_{gen,R+j}-{\overline{x}}_{R+j}\right)+\sum _{j=1,{\tilde{a}}_{k,j}<0}^{N-R}{\tilde{a}}_{k,j}\left(x_{gen,R+j}-{\underline{x}}_{R+j}\right)}{\sum _{j=1,{\tilde{a}}_{i,j}>0}^{N-R}{\tilde{a}}_{i,j}\left(x_{gen,R+j}-{\overline{x}}_{R+j}\right)+\sum _{j=1,{\tilde{a}}_{i,j}<0}^{N-R}{\tilde{a}}_{i,j}\left(x_{gen,R+j}-{\underline{x}}_{R+j}\right)}\times \left[\left(b_{RREF,i}-{\overline{x}}_i\right)-\left(\sum _{j=1,{\tilde{a}}_{k,j}>0}^{N-R}{\tilde{a}}_{i,j}{\overline{x}}_{R+j}+\sum _{j=1,{\tilde{a}}_{k,j}<0}^{N-R}{\tilde{a}}_{i,j}{\underline{x}}_{R+j}+\sum _{j=1,{\tilde{a}}_{k,j}=0}^{N-R}{\tilde{a}}_{i,j}{x}_{gen,R+j}\right)\right]\\ {\overline{\tilde{b}}}_{cor,k}^{\prime \prime }\left(i\right)& =\left(\sum _{j=1,{\tilde{a}}_{k,j}>0}^{N-R}{\tilde{a}}_{i,j}{\overline{x}}_{R+j}+\sum _{j=1,{\tilde{a}}_{k,j}<0}^{N-R}{\tilde{a}}_{i,j}{\underline{x}}_{R+j}\right)+\frac{\sum _{j=1,{\tilde{a}}_{k,j}>0}^{N-R}{\tilde{a}}_{k,j}\left(x_{gen,R+j}-{\overline{x}}_{R+j}\right)+\sum _{j=1,{\tilde{a}}_{k,j}<0}^{N-R}{\tilde{a}}_{k,j}\left(x_{gen,R+j}-{\underline{x}}_{R+j}\right)}{\sum _{j=1,{\tilde{a}}_{i,j}>0}^{N-R}{\tilde{a}}_{i,j}\left(x_{gen,R+j}-{\overline{x}}_{R+j}\right)+\sum _{j=1,{\tilde{a}}_{i,j}<0}^{N-R}{\tilde{a}}_{i,j}\left(x_{gen,R+j}-{\underline{x}}_{R+j}\right)}\times \left[\left(b_{RREF,i}-{\underline{x}}_i\right)-\left(\sum _{j=1,{\tilde{a}}_{k,j}>0}^{N-R}{\tilde{a}}_{i,j}{\overline{x}}_{R+j}+\sum _{j=1,{\tilde{a}}_{k,j}<0}^{N-R}{\tilde{a}}_{i,j}{\underline{x}}_{R+j}+\sum _{j=1,{\tilde{a}}_{k,j}=0}^{N-R}{\tilde{a}}_{i,j}{x}_{gen,R+j}\right)\right]\end{array}\end{array}\text{\hspace{1em}(28)}$$

    情况2.2.2.2.2.2.2：当$\sum _{j=1,{\tilde{a}}_{k,j}>0}^{N-R}{\tilde{a}}_{k,j}\left(x_{gen,R+j}-{\overline{x}}_{R+j}\right)+\sum _{j=1,{\tilde{a}}_{k,j}<0}^{N-R}{\tilde{a}}_{k,j}\left(x_{gen,R+j}-{\underline{x}}_{R+j}\right)$与$\sum _{j=1,{\tilde{a}}_{i,j}>0}^{N-R}{\tilde{a}}_{i,j}\left(x_{gen,R+j}-{\overline{x}}_{R+j}\right)+\sum _{j=1,{\tilde{a}}_{i,j}<0}^{N-R}{\tilde{a}}_{i,j}\left(x_{gen,R+j}-{\underline{x}}_{R+j}\right)$异号时，${\underline{\tilde{b}}}_{cor,k}^{\prime }\left(i\right)$、${\overline{\tilde{b}}}_{cor,k}^{\prime }\left(i\right)$、${\underline{\tilde{b}}}_{cor,k}^{\prime \prime }\left(i\right)$、${\overline{\tilde{b}}}_{cor,k}^{\prime \prime }\left(i\right)$分别满足
$$\begin{array}{c}\begin{array}{rl}{\underline{\tilde{b}}}_{cor,k}^{\prime }\left(i\right)& =\left(\sum _{j=1,{\tilde{a}}_{k,j}>0}^{N-R}{\tilde{a}}_{i,j}{\overline{x}}_{R+j}+\sum _{j=1,{\tilde{a}}_{k,j}<0}^{N-R}{\tilde{a}}_{i,j}{\underline{x}}_{R+j}\right)+\frac{\sum _{j=1,{\tilde{a}}_{k,j}>0}^{N-R}{\tilde{a}}_{k,j}\left(x_{gen,R+j}-{\underline{x}}_{R+j}\right)+\sum _{j=1,{\tilde{a}}_{k,j}<0}^{N-R}{\tilde{a}}_{k,j}\left(x_{gen,R+j}-{\overline{x}}_{R+j}\right)}{\sum _{j=1,{\tilde{a}}_{k,j}>0}^{N-R}{\tilde{a}}_{i,j}\left(x_{gen,R+j}-{\underline{x}}_{R+j}\right)+\sum _{j=1,{\tilde{a}}_{k,j}<0}^{N-R}{\tilde{a}}_{i,j}\left(x_{gen,R+j}-{\overline{x}}_{R+j}\right)}\times \left[\left(\sum _{j=1,{\tilde{a}}_{i,j}>0}^{N-R}{\tilde{a}}_{i,j}{\overline{x}}_{R+j}+\sum _{j=1,{\tilde{a}}_{i,j}<0}^{N-R}{\tilde{a}}_{i,j}{\underline{x}}_{R+j}\right)-\left(\sum _{j=1,{\tilde{a}}_{k,j}>0}^{N-R}{\tilde{a}}_{i,j}{\underline{x}}_{R+j}+\sum _{j=1,{\tilde{a}}_{k,j}<0}^{N-R}{\tilde{a}}_{i,j}{\overline{x}}_{R+j}+\sum _{j=1,{\tilde{a}}_{k,j}=0}^{N-R}{\tilde{a}}_{i,j}{x}_{gen,R+j}\right)\right]\\ {\overline{\tilde{b}}}_{cor,k}^{\prime }\left(i\right)& =\left(\sum _{j=1,{\tilde{a}}_{k,j}>0}^{N-R}{\tilde{a}}_{i,j}{\overline{x}}_{R+j}+\sum _{j=1,{\tilde{a}}_{k,j}<0}^{N-R}{\tilde{a}}_{i,j}{\underline{x}}_{R+j}\right)+\frac{\sum _{j=1,{\tilde{a}}_{k,j}>0}^{N-R}{\tilde{a}}_{k,j}\left(x_{gen,R+j}-{\underline{x}}_{R+j}\right)+\sum _{j=1,{\tilde{a}}_{k,j}<0}^{N-R}{\tilde{a}}_{k,j}\left(x_{gen,R+j}-{\overline{x}}_{R+j}\right)}{\sum _{j=1,{\tilde{a}}_{k,j}>0}^{N-R}{\tilde{a}}_{i,j}\left(x_{gen,R+j}-{\underline{x}}_{R+j}\right)+\sum _{j=1,{\tilde{a}}_{k,j}<0}^{N-R}{\tilde{a}}_{i,j}\left(x_{gen,R+j}-{\overline{x}}_{R+j}\right)}\times \left[\left(\sum _{j=1,{\tilde{a}}_{i,j}>0}^{N-R}{\tilde{a}}_{i,j}{\underline{x}}_{R+j}+\sum _{j=1,{\tilde{a}}_{i,j}<0}^{N-R}{\tilde{a}}_{i,j}{\overline{x}}_{R+j}\right)-\left(\sum _{j=1,{\tilde{a}}_{k,j}>0}^{N-R}{\tilde{a}}_{i,j}{\underline{x}}_{R+j}+\sum _{j=1,{\tilde{a}}_{k,j}<0}^{N-R}{\tilde{a}}_{i,j}{\overline{x}}_{R+j}+\sum _{j=1,{\tilde{a}}_{k,j}=0}^{N-R}{\tilde{a}}_{i,j}{x}_{gen,R+j}\right)\right]\\ {\underline{\tilde{b}}}_{cor,k}^{\prime \prime }\left(i\right)& =\left(\sum _{j=1,{\tilde{a}}_{k,j}>0}^{N-R}{\tilde{a}}_{i,j}{\overline{x}}_{R+j}+\sum _{j=1,{\tilde{a}}_{k,j}<0}^{N-R}{\tilde{a}}_{i,j}{\underline{x}}_{R+j}\right)+\frac{\sum _{j=1,{\tilde{a}}_{k,j}>0}^{N-R}{\tilde{a}}_{k,j}\left(x_{gen,R+j}-{\overline{x}}_{R+j}\right)+\sum _{j=1,{\tilde{a}}_{k,j}<0}^{N-R}{\tilde{a}}_{k,j}\left(x_{gen,R+j}-{\underline{x}}_{R+j}\right)}{\sum _{j=1,{\tilde{a}}_{i,j}>0}^{N-R}{\tilde{a}}_{i,j}\left(x_{gen,R+j}-{\overline{x}}_{R+j}\right)+\sum _{j=1,{\tilde{a}}_{i,j}<0}^{N-R}{\tilde{a}}_{i,j}\left(x_{gen,R+j}-{\underline{x}}_{R+j}\right)}\times \left[\left(b_{RREF,i}-{\underline{x}}_i\right)-\left(\sum _{j=1,{\tilde{a}}_{k,j}>0}^{N-R}{\tilde{a}}_{i,j}{\overline{x}}_{R+j}+\sum _{j=1,{\tilde{a}}_{k,j}<0}^{N-R}{\tilde{a}}_{i,j}{\underline{x}}_{R+j}+\sum _{j=1,{\tilde{a}}_{k,j}=0}^{N-R}{\tilde{a}}_{i,j}{x}_{gen,R+j}\right)\right]\\ {\overline{\tilde{b}}}_{cor,k}^{\prime \prime }\left(i\right)& =\left(\sum _{j=1,{\tilde{a}}_{k,j}>0}^{N-R}{\tilde{a}}_{i,j}{\overline{x}}_{R+j}+\sum _{j=1,{\tilde{a}}_{k,j}<0}^{N-R}{\tilde{a}}_{i,j}{\underline{x}}_{R+j}\right)+\frac{\sum _{j=1,{\tilde{a}}_{k,j}>0}^{N-R}{\tilde{a}}_{k,j}\left(x_{gen,R+j}-{\overline{x}}_{R+j}\right)+\sum _{j=1,{\tilde{a}}_{k,j}<0}^{N-R}{\tilde{a}}_{k,j}\left(x_{gen,R+j}-{\underline{x}}_{R+j}\right)}{\sum _{j=1,{\tilde{a}}_{i,j}>0}^{N-R}{\tilde{a}}_{i,j}\left(x_{gen,R+j}-{\overline{x}}_{R+j}\right)+\sum _{j=1,{\tilde{a}}_{i,j}<0}^{N-R}{\tilde{a}}_{i,j}\left(x_{gen,R+j}-{\underline{x}}_{R+j}\right)}\times \left[\left(b_{RREF,i}-{\overline{x}}_i\right)-\left(\sum _{j=1,{\tilde{a}}_{k,j}>0}^{N-R}{\tilde{a}}_{i,j}{\overline{x}}_{R+j}+\sum _{j=1,{\tilde{a}}_{k,j}<0}^{N-R}{\tilde{a}}_{i,j}{\underline{x}}_{R+j}+\sum _{j=1,{\tilde{a}}_{k,j}=0}^{N-R}{\tilde{a}}_{i,j}{x}_{gen,R+j}\right)\right]\end{array}\end{array}\text{\hspace{1em}(29)}$$

    将$x_{cor,R+1},\cdots ,x_{cor,N}$代入式(13)，则修正后的$x_{cor,1},\cdots ,x_{cor,R}$满足
x c o r , i = b R R E F , i − ∑ j = 1 N − R a ~ i , j x c o r , R + j , i ∈ { 1 , ⋯   , R } \begin{equation}\tag{30} x_{cor, i} = b_{RREF, i} - \displaystyle \sum_{j = 1}^{N - R} \tilde{a}_{i, j} x_{cor,R + j}, i ∈ \{ 1, \cdots, R \} \end{equation} xcor,i​=bRREF,i​−j=1∑N−R​a~i,j​xcor,R+j​,i∈{1,⋯,R}​(30)​
    显然，若${\tilde{b}}_{cor,k}$不满足式(25)所示的约束条件，则必存在 i ∈ { 1 , ⋯   , R } i ∈ \{ 1, \cdots, R \} i∈{1,⋯,R}，使得$x_{cor,i}$不能满足式(1)中的上下界约束，进而无法将${\mathbf{x}}_{\mathbf{gen}}$修正为满足式(1)中所有约束条件的${\mathbf{x}}_{\mathbf{cor}}$。

MATLAB代码

linear_equations_constraints_correction函数

function [x_cor, correction_info] = linear_equations_constraints_correction(x_gen, lb, ub, A_eq, b_eq, ...
                                    max_runtime)
    % LINEAR_EQUATIONS_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('max_runtime', 'var')
        max_runtime = ceil(1000 / size(x_gen, 2));
    end
    
    Ab_eq = [A_eq, b_eq];
    
    M = size(A_eq, 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 = sum(any(Ab_eq_RREF(:, 1:end - 1), 2), 1);
    R_Ab = sum(any(Ab_eq_RREF, 2), 1);
    
    if R_A < R_Ab
        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;
    
    else
        if R_A == N_D
            x_cor = Ab_eq_RREF(1:N_D, end);
    
        else
            x_cor = x_gen;

            if R_A == 1
                A_eq_candi = A_eq(1, :);
                Eq_gen = A_eq_candi * x_gen;

                b_eq_candi_lb = A_eq_candi * ((A_eq_candi > 0)' .* lb + (A_eq_candi < 0)' .* ub);
                b_eq_candi_ub = A_eq_candi * ((A_eq_candi > 0)' .* ub + (A_eq_candi < 0)' .* lb);
                delta_plus = A_eq_candi * ((A_eq_candi > 0)' .* (x_gen - lb) + (A_eq_candi < 0)' .* (x_gen - ub));
                delta_minus = A_eq_candi * ((A_eq_candi > 0)' .* (x_gen - ub) + (A_eq_candi < 0)' .* (x_gen - lb));

                for indi = 1:N_P
                    if Eq_gen(:, indi) > b_eq
                        x_cor(A_eq_candi > 0, indi) = (((b_eq - b_eq_candi_lb) / delta_plus(indi)) * (x_gen(A_eq_candi > 0, indi) - lb(A_eq_candi > 0)) ...
                                                      + lb(A_eq_candi > 0))';
                        x_cor(A_eq_candi == 0, indi) = (x_gen(A_eq_candi == 0, indi))';
                        x_cor(A_eq_candi < 0, indi) = (((b_eq - b_eq_candi_lb) / delta_plus(indi)) * (x_gen(A_eq_candi < 0, indi) - ub(A_eq_candi < 0)) ...
                                                      + ub(A_eq_candi < 0))';
                    elseif Eq_gen(:, indi) < b_eq
                        x_cor(A_eq_candi > 0, indi) = (((b_eq - b_eq_candi_ub) / delta_minus(indi)) * (x_gen(A_eq_candi > 0, indi) - ub(A_eq_candi > 0)) ...
                                                      + ub(A_eq_candi > 0))';
                        x_cor(A_eq_candi == 0, indi) = (x_gen(A_eq_candi == 0, indi))';
                        x_cor(A_eq_candi < 0, indi) = (((b_eq - b_eq_candi_ub) / delta_minus(indi)) * (x_gen(A_eq_candi < 0, indi) - lb(A_eq_candi < 0)) ...
                                                      + lb(A_eq_candi < 0))';
                    end
                end

            else
                lb_base = lb(R_A + 1:end, :);
                ub_base = ub(R_A + 1:end, :);
                A_eq_tilde = Ab_eq_RREF(1:R_A, R_A + 1:end - 1);
                b_eq_tilde_lb = diag(A_eq_tilde * ((A_eq_tilde > 0)' .* lb_base + (A_eq_tilde < 0)' .* ub_base));
                b_eq_tilde_ub = diag(A_eq_tilde * ((A_eq_tilde > 0)' .* ub_base + (A_eq_tilde < 0)' .* lb_base));
                b_eq_tilde_min = Ab_eq_RREF(:, end) - ub(1:R_A, end);
                b_eq_tilde_max = Ab_eq_RREF(:, end) - lb(1:R_A, end);
                
                if any(b_eq_tilde_max < b_eq_tilde_lb) || any(b_eq_tilde_min > b_eq_tilde_ub)
                    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;
                else
                    run = 0;
                    while 1
                        x_cor = x_gen;
                        x_gen_base = x_gen(R_A + 1:end, :);
                        Eq_tilde_gen = A_eq_tilde * x_gen_base;
                        does_A_eq_tilde_x_gen_base_exceed = (Eq_tilde_gen < b_eq_tilde_min) ...
                                                            | (Eq_tilde_gen > b_eq_tilde_max);

                        for indi = 1:N_P
                            if any(does_A_eq_tilde_x_gen_base_exceed(:, indi))
                                flag_A_eq_tilde_x_gen_base_exceed = find(does_A_eq_tilde_x_gen_base_exceed(:, indi));
                                row_A_eq_tilde_candi = flag_A_eq_tilde_x_gen_base_exceed(randi(length(flag_A_eq_tilde_x_gen_base_exceed)));
                                A_eq_tilde_candi = A_eq_tilde(row_A_eq_tilde_candi, :);
                                Eq_tilde_candi_gen = A_eq_tilde_candi * x_gen_base(:, indi);
                                b_eq_k_tilde_lb = A_eq_tilde_candi * ((A_eq_tilde_candi > 0)' .* lb_base ...
                                                  + (A_eq_tilde_candi < 0)' .* ub_base);
                                b_eq_k_tilde_ub = A_eq_tilde_candi * ((A_eq_tilde_candi > 0)' .* ub_base ...
                                                  + (A_eq_tilde_candi < 0)' .* lb_base);
                                delta_plus = A_eq_tilde_candi * ((A_eq_tilde_candi > 0)' .* (x_gen_base(:, indi) - lb_base) ...
                                             + (A_eq_tilde_candi < 0)' .* (x_gen_base(:, indi) - ub_base));
                                delta_minus = A_eq_tilde_candi * ((A_eq_tilde_candi > 0)' .* (x_gen_base(:, indi) - ub_base) ...
                                             + (A_eq_tilde_candi < 0)' .* (x_gen_base(:, indi) - lb_base));
                                x_cor_base = zeros(N_D - R_A, 1);
        
                                if Eq_tilde_candi_gen > b_eq_tilde_max(row_A_eq_tilde_candi)
                                    b_eq_k_tilde_best = [b_eq_k_tilde_lb + (delta_plus * (b_eq_tilde_min - A_eq_tilde * ((A_eq_tilde_candi > 0)' ...
                                                        .* lb_base + (A_eq_tilde_candi < 0)' .* ub_base + (A_eq_tilde_candi == 0)' .* x_gen_base(:, indi)))) ...
                                                        ./ (A_eq_tilde * ((A_eq_tilde_candi > 0)' .* (x_gen_base(:, indi) - lb_base) ...
                                                        + (A_eq_tilde_candi < 0)' .* (x_gen_base(:, indi) - ub_base))), ...
                                                        b_eq_k_tilde_lb + (delta_plus * (b_eq_tilde_max - A_eq_tilde * ((A_eq_tilde_candi > 0)' ...
                                                        .* lb_base + (A_eq_tilde_candi < 0)' .* ub_base + (A_eq_tilde_candi == 0)' .* x_gen_base(:, indi)))) ...
                                                        ./ (A_eq_tilde * ((A_eq_tilde_candi > 0)' .* (x_gen_base(:, indi) - lb_base) ...
                                                        + (A_eq_tilde_candi < 0)' .* (x_gen_base(:, indi) - ub_base)));
                                                        b_eq_k_tilde_lb, b_eq_k_tilde_ub];
                                    b_eq_k_tilde_min = max(min(b_eq_k_tilde_best, [], 2), [], 1);
                                    b_eq_k_tilde_max = min(max(b_eq_k_tilde_best, [], 2), [], 1);
                                    if b_eq_k_tilde_min <= b_eq_k_tilde_max
                                        b_eq_k_tilde = b_eq_k_tilde_min + rand * (b_eq_k_tilde_max - b_eq_k_tilde_min);
                                    else
                                        b_eq_k_tilde = b_eq_k_tilde_best(row_A_eq_tilde_candi, 1) + rand * (b_eq_k_tilde_best(row_A_eq_tilde_candi, 2)...
                                                       - b_eq_k_tilde_best(row_A_eq_tilde_candi, 1));
                                    end
                                    x_cor_base(A_eq_tilde_candi > 0) = (((b_eq_k_tilde - b_eq_k_tilde_lb) / delta_plus) ...
                                                                       * (x_gen_base(A_eq_tilde_candi > 0, indi) - lb_base(A_eq_tilde_candi > 0)) ...
                                                                       + lb_base(A_eq_tilde_candi > 0))';
                                    x_cor_base(A_eq_tilde_candi == 0) = (x_gen_base(A_eq_tilde_candi == 0))';
                                    x_cor_base(A_eq_tilde_candi < 0) = (((b_eq_k_tilde - b_eq_k_tilde_lb) / delta_plus) ...
                                                                       * (x_gen_base(A_eq_tilde_candi < 0, indi) - ub_base(A_eq_tilde_candi < 0)) ...
                                                                       + ub_base(A_eq_tilde_candi < 0))';
                                elseif Eq_tilde_candi_gen < b_eq_tilde_min(row_A_eq_tilde_candi)
                                    b_eq_k_tilde_best = [b_eq_k_tilde_ub + (delta_minus * (b_eq_tilde_min - A_eq_tilde * ((A_eq_tilde_candi > 0)' .* ub_base ...
                                                        + (A_eq_tilde_candi < 0)' .* lb_base + (A_eq_tilde_candi == 0)' .* x_gen_base(:, indi))) ./ (A_eq_tilde ...
                                                        * ((A_eq_tilde_candi > 0)' .* (x_gen_base(:, indi) - ub_base) + (A_eq_tilde_candi < 0)' .* ...
                                                        (x_gen_base(:, indi) - lb_base)))), ...
                                                        b_eq_k_tilde_ub + (delta_minus * (b_eq_tilde_max - A_eq_tilde * ((A_eq_tilde_candi > 0)' .* ub_base ...
                                                        + (A_eq_tilde_candi < 0)' .* lb_base + (A_eq_tilde_candi == 0)' .* x_gen_base(:, indi))) ./ (A_eq_tilde ...
                                                        * ((A_eq_tilde_candi > 0)' .* (x_gen_base(:, indi) - ub_base) + (A_eq_tilde_candi < 0)' .* ...
                                                        (x_gen_base(:, indi) - lb_base))));
                                                        b_eq_k_tilde_lb, b_eq_k_tilde_ub];
                                    b_eq_k_tilde_min = max(min(b_eq_k_tilde_best, [], 2), [], 1);
                                    b_eq_k_tilde_max = min(max(b_eq_k_tilde_best, [], 2), [], 1);
                                    if b_eq_k_tilde_min <= b_eq_k_tilde_max
                                        b_eq_k_tilde = b_eq_k_tilde_min + rand * (b_eq_k_tilde_max - b_eq_k_tilde_min);
                                    else
                                        b_eq_k_tilde = b_eq_k_tilde_best(row_A_eq_tilde_candi, 1) + rand * (b_eq_k_tilde_best(row_A_eq_tilde_candi, 2)...
                                                       - b_eq_k_tilde_best(row_A_eq_tilde_candi, 1));
                                    end
                                    x_cor_base(A_eq_tilde_candi > 0) = (((b_eq_k_tilde - b_eq_k_tilde_ub) / delta_minus) ...
                                                                       * (x_gen_base(A_eq_tilde_candi > 0, indi) - ub_base(A_eq_tilde_candi > 0)) ...
                                                                       + ub_base(A_eq_tilde_candi > 0))';
                                    x_cor_base(A_eq_tilde_candi == 0) = (x_gen_base(A_eq_tilde_candi == 0))';
                                    x_cor_base(A_eq_tilde_candi < 0) = (((b_eq_k_tilde - b_eq_k_tilde_ub) / delta_minus) ...
                                                                       * (x_gen_base(A_eq_tilde_candi < 0, indi) - lb_base(A_eq_tilde_candi < 0)) ...
                                                                       + lb_base(A_eq_tilde_candi < 0))';
                                end
                                x_cor(R_A + 1:end, indi) = x_cor_base;
                            end
                        end
                        x_cor(1:R_A, :) = Ab_eq_RREF(1:R_A, end) - A_eq_tilde * x_cor(R_A + 1: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(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);

                        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            

        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_equations_constraints_correction(x_gen, lb, ub, A_eq, b_eq);

Problem 2

    对应博客 《一种多线性方程约束下的生成随机数修正方法（上）》实例分析部分中的实例1。

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_equations_constraints_correction(x_gen, lb, ub, A_eq, b_eq);

Problem 3

    对应博客 《一种多线性方程约束下的生成随机数修正方法（上）》实例分析部分中的实例2。

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; -7, -4, 8, -5, 0, 1];
b_eq = [-2; 14];

[x_cor, correction_info] = linear_equations_constraints_correction(x_gen, lb, ub, A_eq, b_eq);

Problem 4

    对应博客 《一种多线性方程约束下的生成随机数修正方法（上）》和博客 《一种多线性方程约束下的生成随机数修正方法（下）》实例分析部分中的实例3。

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; 10, 3, -3, 6, -7, 2; 2, -3, -7, 5, -6, 3];
b_eq = [-2; 14; 9; -8];

[x_cor, correction_info] = linear_equations_constraints_correction(x_gen, lb, ub, A_eq, b_eq);

研究目标

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