1. Optimization and linear programming
The three elements of optimization problem are decision variables, objective function and constraints.
Linear programming is an optimization method to study the extreme value problem of linear objective function under linear constraints. It is often used to solve the problem of using existing resources to obtain the optimal decision.
Simple linear programming problems can be solved by Lingo software. Matlab and Python also have library functions or solvers for solving linear programming problems, which are easy to learn and use, and do not need simulated annealing algorithm. However, integer programming, mixed programming, 0 / 1 programming, quadratic programming, nonlinear programming and combinatorial optimization problems derived from general linear programming problems can not be handled by calling a library function. The simulated annealing algorithm has good adaptability in many complex problems, and can be used as an entrylevel general intelligent algorithm to learn.
In other words, if we only deal with linear programming problems, we should not use simulated annealing algorithm. However, if it is a complex optimization problem that cannot be handled by the existing methods, or you don't know what method to deal with a certain kind of optimization problem, simulated annealing algorithm can still be used to solve it.
In this paper, the penalty function method is used to analyze the simulated annealing algorithm to deal with linear programming problems. The related contents are also applicable to nonlinear programming problems.
2. Simulated annealing algorithm deals with constraints
Linear programming problem is a constrained optimization problem, and simulated annealing algorithm is more suitable to deal with unconstrained optimization problems. For the constraints in optimization problems, simulated annealing algorithm has several common methods:

1. Upper and lower limits of decision variables.
Such constraints are relatively easy to handle, and can be solved as long as the initial solution and new solution are set between the upper and lower limits of the value of decision variables. For example, (1) set the upper and lower limits of the random number generating the new solution as the upper and lower limits of the decision variable, i.e. [Xmin, Xmax]; (2) Set the upper and lower limits of the random number generating the new solution as the upper and lower limits of the current solution and the decision variable, i.e. [Xnow, Xmax]; (3) Through conditional judgment, when the new solution exceeds the upper and lower limits of the decision variable, make it take the upper and lower limits, that is, xNew = max(min(xNew, xMax), xMin). Of course, these processing methods will affect the probability distribution of random numbers and thus the optimization performance of simulated annealing algorithm, which will not be discussed in depth here. 
2. Test method to deal with inequality constraints.
In the iterative process of simulated annealing algorithm, the new solution generated each time is substituted into each inequality constraint function to judge whether the constraint conditions are met; If the new solution does not meet the constraints, discard the new solution and return to generate a new solution for inspection until the new solution meets all the constraints. The idea of this method is simple, and each iteration is carried out in the feasible region. However, for complex problems with many and harsh constraints, the new solutions generated many times can not meet the constraints, which will make the calculation time very long or even stagnate. 
3. The elimination method deals with equality constraint problems.
For equality constraints, it is difficult to meet the constraints through the randomly generated new solutions, and the test method can not be used generally. The elimination method is to solve the equation, express a decision variable in the equality constraint as the linear relationship of other decision variables, and then substitute it into the objective function and inequality constraint conditions, so as to eliminate the constraint conditions. The elimination method not only solves the equality constraint problem, but also reduces the number of decision variables, which effectively simplifies the complexity of the optimization problem. It is a method to kill two birds with one stone. However, for nonlinear equality constraints, it is also very difficult to solve nonlinear equations, and the elimination method is not generally applicable. 
4. A more general method to deal with constraints is the penalty function method, which is introduced below. YouCans, XUPT
3. Penalty function method
Penalty function method is a kind of commonly used technology to deal with constraints. It is very effective to deal with constraints in simulated annealing algorithm. The idea of the method is to transform the constraint condition into a penalty function and add it to the original objective function to construct a new objective function; When the constraint conditions are not satisfied, the new objective function is discarded by making the penalty function worse.
Penalty function method includes outer point method and inner point method. The outer point method penalizes the points outside the feasible region (i.e. the points that do not meet the constraints) and does not penalize the points inside the feasible region, so that the iterative points approach the feasible region D. The interior point method searches inside the feasible region, and the constraint boundary acts like a fence, so that the objective function cannot pass through, so the search point is limited to the feasible region, so it is only applicable to inequality constraints.
Consider constrained optimization problems:
m
i
n
:
f
(
X
)
min: f(X)
min:f(X)
Meet the limit
c
i
(
X
)
≤
0
c_i(X)\leq0
ci(X)≤0
The penalty function method transforms the problem into the following unconstrained problem
m
i
n
:
L
k
(
X
)
=
f
(
X
)
+
σ
k
∑
i
g
(
c
i
(
X
)
)
min: L_k(X)=f(X)+\sigma_k\sum_{i}{g(c_i(X))}
min:Lk(X)=f(X)+σki∑g(ci(X))
among
g
(
c
i
(
X
)
)
=
m
a
x
(
0
,
c
i
(
X
)
)
2
g(c_i(X))=max(0,c_i(X))^2
g(ci(X))=max(0,ci(X))2
In the above equation, g ( c i ( X ) ) g(ci(X)) g(ci(X)) is called the external penalty function, s i g m a ( k ) sigma(k) sigma(k) is called penalty factor. In each iteration, s i g m a ( k ) sigma(k) sigma(k) increases, and then the unconstrained problem is solved. The results of each iteration will form a sequence, and the limit of this sequence is the solution of the original constraint problem.
4. Digital analog case
Although simulated annealing algorithm is not recommended for linear programming problems. However, in order to facilitate understanding, this paper still uses the previous linear programming problem as a case to deal with constraints. For nonlinear programming problems and nonlinear constraint problems, the treatment methods are similar, which will be introduced later.
4.1 Problem Description:
A factory produces two kinds of drinks, 6 kg of raw materials and 10 workers per 100 boxes of drinks, making a profit of 100000 yuan; Each hundred boxes of B drinks need 5kg of raw materials and 20 workers, making a profit of 90000 yuan.
Today, the factory has a total of 60kg of raw materials and 150 workers. Due to other conditions, the output of a beverage is no more than 800 boxes.
(1) how to arrange the production plan, that is, how much are the two drinks produced to maximize the profit?
4.2 problem modeling:
Decision variables:
x1: output of a beverage (unit: 100 cases)
x2: output of B beverage (unit: 100 cases)
Objective function:
max fx = 10x1 + 9x2
Constraints:
6x1 + 5x2 <= 60
10x1 + 20x2 <= 150
Value range:
Given conditions: x1, X2 > = 0, X1 < = 8
Derivation conditions: from x1, X2 > = 0 and 10 * X1 + 20 * x2 < = 150, we can know: 0 < = X1 < = 15; 0<=x2<=7.5
Therefore, 0 < = X1 < = 8, 0 < = x2 < = 7.5
4.3 penalty function method for solving constrained optimization problems:
Construct penalty function:
p1 = (max(0, 6*x1+5*x260))**2
p2 = (max(0, 10*x1+20*x2150))**2
Note: if there are equality constraints, such as x1 + 2*x2 = m, it can also be transformed into penalty function:
p3 = (x1+2*x2m)**2
P(x) = p1 + p2 + ...
Construct augmented objective function:
L(x,m(k)) = min(fx) + m(k)*P(x)
m(k): penalty factor, which increases gradually with the number of iterations K
In the simulated annealing algorithm, m(k) increases with the number of iterations of the outer loop, but it should remain unchanged in the inner loop.
5. Simulated annealing algorithm Python program: penalty function method to solve constrained optimization linear programming problem
# Simulated annealing algorithm program: penalty function method to solve linear programming problems # Program: SimulatedAnnealing_v2.py # Purpose: Simulated annealing algorithm for function optimization # v2.0: use the penalty function method to deal with the constraint problem # Copyright 2021 YouCans, XUPT # Crated: 20210501 # * coding: utf8 * import math # Import module import random # Import module import pandas as pd # Import module YouCans, XUPT import numpy as np # Import the module numpy, which is abbreviated as np import matplotlib.pyplot as plt from datetime import datetime # Subroutine: defines the objective function of the optimization problem def cal_Energy(X, nVar, mk): # m(k): penalty factor, which increases gradually with the number of iterations K p1 = (max(0, 6*X[0]+5*X[1]60))**2 p2 = (max(0, 10*X[0]+20*X[1]150))**2 fx = (10*X[0]+9*X[1]) return fx+mk*(p1+p2) # Subroutine: parameter setting of simulated annealing algorithm def ParameterSetting(): cName = "funcOpt" # Define the problem name YouCans, XUPT nVar = 2 # Given the number of arguments, y = f (x1,... Xn) xMin = [0, 0] # Lower limit of given search space, x1_min,..xn_min xMax = [8, 7.5] # Upper limit of given search space, x1_max,..xn_max tInitial = 100.0 # Set initial annealing temperature tFinal = 1 # Set stop temperature alfa = 0.98 # Set cooling parameters, T(k)=alfa*T(k1) meanMarkov = 100 # Markov chain length, that is, the number of inner loop operations scale = 0.5 # Define the search step size, which can be set to a fixed value or gradually reduced return cName, nVar, xMin, xMax, tInitial, tFinal, alfa, meanMarkov, scale # Simulated annealing algorithm def OptimizationSSA(nVar,xMin,xMax,tInitial,tFinal,alfa,meanMarkov,scale): # ======Initialize random number generator====== randseed = random.randint(1, 100) random.seed(randseed) # The random number generator can set the seed or set it to a specified integer # ======Initial solution of stochastic optimization problem====== xInitial = np.zeros((nVar)) # Initialize, create array for v in range(nVar): # random.uniform(min,max) randomly generates a real number in the range of [min,max] xInitial[v] = random.uniform(xMin[v], xMax[v]) # Call sub function cal_Energy calculates the objective function value of the current solution fxInitial = cal_Energy(xInitial, nVar, 1) # m(k): penalty factor, the initial value is 1 # ======Simulated annealing algorithm initialization====== xNew = np.zeros((nVar)) # Initialize, create array xNow = np.zeros((nVar)) # Initialize, create array xBest = np.zeros((nVar)) # Initialize, create array xNow[:] = xInitial[:] # Initialize the current solution and set the initial solution as the current solution xBest[:] = xInitial[:] # Initialize the optimal solution and set the current solution as the optimal solution fxNow = fxInitial # Set the objective function of the initial solution to the current value fxBest = fxInitial # Set the objective function of the current solution as the optimal value print('x_Initial:{:.6f},{:.6f},\tf(x_Initial):{:.6f}'.format(xInitial[0], xInitial[1], fxInitial)) recordIter = [] # Initialization, number of external cycles recordFxNow = [] # Initialization, the objective function value of the current solution recordFxBest = [] # Initialization, the objective function value of the optimal solution recordPBad = [] # Initialization, acceptance probability of inferior solution kIter = 0 # Number of external loop iterations, number of temperature states totalMar = 0 # Total Markov chain length totalImprove = 0 # fxBest improvement times nMarkov = meanMarkov # Fixed length Markov chain # ======Start simulated annealing optimization====== # The external circulation ends when the current temperature reaches the termination temperature tNow = tInitial # Initialize current temperature while tNow >= tFinal: # The external circulation ends when the current temperature reaches the termination temperature # At the current temperature, a sufficient number of state transitions (nMarkov) are performed to achieve thermal equilibrium kBetter = 0 # Number of times to obtain highquality solutions kBadAccept = 0 # Number of times to accept inferior solutions kBadRefuse = 0 # Number of rejections of inferior solutions # Inner loop, the number of cycles is the length of Markov chain for k in range(nMarkov): # Inner loop, the number of cycles is the length of Markov chain totalMar += 1 # Total Markov chain length counter # Generate new solutions # Generate a new solution: generate a new solution by randomly disturbing near the current solution. The new solution must be within the range of [min,max] # Scheme 1: only one of the n variables is disturbed, and the other n1 variables remain unchanged xNew[:] = xNow[:] v = random.randint(0, nVar1) # Generate random numbers between [0,nVar1] xNew[v] = xNow[v] + scale * (xMax[v]xMin[v]) * random.normalvariate(0, 1) # random. Normal variable (0, 1): generate a random real number that obeys a normal distribution with a mean of 0 and a standard deviation of 1 xNew[v] = max(min(xNew[v], xMax[v]), xMin[v]) # Ensure that the new solution is within the range of [min,max] # Calculate the objective function and energy difference # Call sub function cal_Energy calculates the objective function value of the new solution fxNew = cal_Energy(xNew, nVar, kIter) deltaE = fxNew  fxNow # Accept the new solution according to Metropolis guidelines # Acceptance judgment: decide whether to accept the new solution according to Metropolis criteria if fxNew < fxNow: # Better solution: if the objective function of the new solution is better than the current solution, the new solution is accepted accept = True kBetter += 1 else: # Tolerant solution: if the objective function of the new solution is worse than the current solution, the new solution will be accepted with a certain probability pAccept = math.exp(deltaE / tNow) # Calculate the state transition probability of the tolerant solution if pAccept > random.random(): accept = True # Accept inferior solutions kBadAccept += 1 else: accept = False # Reject inferior solutions kBadRefuse += 1 # Save new solution if accept == True: # If the new solution is accepted, the new solution is saved as the current solution xNow[:] = xNew[:] fxNow = fxNew if fxNew < fxBest: # If the objective function of the new solution is better than the optimal solution, the new solution is saved as the optimal solution fxBest = fxNew xBest[:] = xNew[:] totalImprove += 1 scale = scale*0.99 # Variable search step size, gradually reduce the search range and improve the search accuracy # Data sorting after internal circulation # Complete the search of current temperature, save data and output pBadAccept = kBadAccept / (kBadAccept + kBadRefuse) # Acceptance probability of inferior solution recordIter.append(kIter) # Current external circulation times recordFxNow.append(round(fxNow, 4)) # Objective function value of current solution recordFxBest.append(round(fxBest, 4)) # Objective function value of optimal solution recordPBad.append(round(pBadAccept, 4)) # Objective function value of optimal solution if kIter%10 == 0: # Modular operation, remainder of quotient print('i:{},t(i):{:.2f}, badAccept:{:.6f}, f(x)_best:{:.6f}'.\ format(kIter, tNow, pBadAccept, fxBest)) # Slowly cool down to a new temperature. Cooling curve: T(k)=alfa*T(k1) tNow = tNow * alfa kIter = kIter + 1 fxBest = cal_Energy(xBest, nVar, kIter) # Since the penalty factor increases after iteration, the augmented objective function needs to be reconstructed # ======End the simulated annealing process====== print('improve:{:d}'.format(totalImprove)) return kIter,xBest,fxBest,fxNow,recordIter,recordFxNow,recordFxBest,recordPBad # Result verification and output def ResultOutput(cName,nVar,xBest,fxBest,kIter,recordFxNow,recordFxBest,recordPBad,recordIter): # ======Verification and output of optimization results====== fxCheck = cal_Energy(xBest, nVar, kIter) if abs(fxBest  fxCheck)>1e3: # Objective function test print("Error 2: Wrong total millage!") return else: print("\nOptimization by simulated annealing algorithm:") for i in range(nVar): print('\tx[{}] = {:.6f}'.format(i,xBest[i])) print('\n\tf(x):{:.6f}'.format(cal_Energy(xBest,nVar,0))) return # main program def main(): # YouCans, XUPT # Parameter setting, optimization problem parameter definition, simulated annealing algorithm parameter setting [cName, nVar, xMin, xMax, tInitial, tFinal, alfa, meanMarkov, scale] = ParameterSetting() # print([nVar, xMin, xMax, tInitial, tFinal, alfa, meanMarkov, scale]) # Simulated annealing algorithm [giter, xbest, fxbest, fxnow, recorditer, recordfxnow, recordfxbest, recordpbad]\ = OptimizationSSA(nVar,xMin,xMax,tInitial,tFinal,alfa,meanMarkov,scale) # print(kIter, fxNow, fxBest, pBadAccept) # Result verification and output ResultOutput(cName, nVar,xBest,fxBest,kIter,recordFxNow,recordFxBest,recordPBad,recordIter) if __name__ == '__main__': main()
6. Operation results
Optimization by simulated annealing algorithm: x[0] = 6.577964 x[1] = 4.111469 f(x):102.782857
reference:
(1) Hu Shanying, Chen Bingzhen, simulated annealing method for global optimization of nonlinear programming problems, Journal of Tsinghua University, 1997,37 (6): 59
(2) Li Qiqiang, simulated annealing algorithm with constraint guidance, systems engineering, 2001,19 (3): 4955
Copyright description:
Original works
Copyright 2021 YouCans, XUPT
Crated: 20210501