static double perturbation(double stest, double dUtempt, double dLtempt, double Para_Perturb) {
        double k1 = 0;
        double k2 = 0;
        double k5;

        // Method 1:
        // Standard Gaussian random number based upon Numerical recipes gasdev and
        // Marsaglia-Bray Algorithm
        double k3 = 2.0;
        while (k3 >= 1.0 || k3 == 0.0) {
            k1 = 2.0 * Math.random() - 1.0;
            k2 = 2.0 * Math.random() - 1.0;
            k3 = k1 * k1 + k2 * k2;
        }
        
        k3 = Math.pow((-2.0 * Math.log(k3) / k3), 0.5);
        double k4 = Math.random();
        if (k4 < 0.5) {
            k5 = k1 * k3;
        } else {
            k5 = k2 * k3;
        }

        // DDS perturbation parameter
        double s_new = stest + Para_Perturb * (dUtempt - dLtempt) * k5;

        // Generate normally distributed random number by MATLAB function
        // s_new = stest + k6 * (dUtempt-dLtempt) * normrnd(0,1) ;
        // Check if s_new is overshooting the feasible bound
        if (s_new > dUtempt) {            	// Upper bound management
            s_new = dUtempt - (s_new - dUtempt);
            if (s_new < dLtempt) {
                s_new = dUtempt;
            }
        } else if (s_new < dLtempt) {      	// Lower bound management
            s_new = dLtempt + (dLtempt - s_new);
            if (s_new > dUtempt) {
                s_new = dLtempt;
            }
        }

        //System.out.println(k5);
        return s_new;
    }