1. Introduction
The development of commercially-off-the-shelf (COTS) components has led to the capability of developing compact spacecraft platforms by combining subsystems from the readily available COTS inventory. For these COTS assembled spacecrafts, the availability of rideshare opportunities on small spacecraft launchers and the prospect of being launched as a secondary payload on larger launch vehicles has led to the rapid advancement of COTS components and expansion of the small spacecraft market. As such, the buses for these small spacecrafts are designed as modular platforms and the components designed to be as diverse and capable as possible, which can be assembled with a payload in a short amount of time, ready to fly [
1]. One such platform is the CubeSat platform, which weighs only a few kilograms and is based on a form factor of a 10 × 10 × 10 cm
3. CubeSats can be composed of a single cube (a “1U” CubeSat) or several cubes combined, forming, for instance, 3U or 6U units [
2]. Although the COTS components are readily available, designing these spacecrafts is a long, expensive endeavor, where a team of experts work together to design the system for a specific mission using their engineering judgement. An alternative design approach is required that shortens and simplifies the spacecraft design process. In this paper, we propose a tool for automated design of a spacecraft using machine learning algorithms that utilize modular, interchangeable components from a COTS inventory to develop a spacecraft design. In our approach, the spacecraft design is presented in the form of a gene that selects subsystem components from an inventory, and a numerical goal function is specified along with constraints, based on which the spacecraft design is evaluated. An evolutionary algorithm (EA) is used to generate a population of candidate solutions, which undergo evolutionary operations such as selection, crossover, and mutation, and the fittest individuals are carried to the next generation, while the unfit individuals are culled off. This method of evolution is performed over hundreds of generations until the fittest individual in the population meets a desired performance metric [
3]. This approach can be very useful to a design team to make informed system design decisions by finding a set of near-optimal solutions that represent the trade-offs among conflicting objectives, such as mass, cost, and performance.
Automated design of engineering systems is not new. When compared to human created designs, computationally derived evolutionary designs have shown competitive advantages in terms of robustness, performance, and creativity. Design optimization through evolutionary computation has been studied for decades. Work has been done to design several satellite subsystems. One such example is an X-band antenna for NASA’s Space Technology 5 (ST5) spacecraft, which is in fact the first computer-evolved hardware in space [
4]. Other work includes design optimization of a power subsystem [
5], low-thrust orbit transfer [
6], spacecraft structure [
7], and a corner retroreflector for satellites [
8], which has shown potential improvements over traditional methods. An evolutionary approach has also been applied to design satellite constellations for continuous regional coverage [
9], nontraditional constellations such as heterogeneous and rideshare constellations using COTS components to increase resiliency and responsiveness [
10], reconfigurable satellite constellations [
11], regional navigation satellite system constellations [
12], pattern forming tasks using robot teams [
13] and even teams of excavating robots for lunar base preparation [
14], which have produced human-competitive designs. This process has also been successfully applied to the design of mobile robots [
15] and water desalination systems [
16]. Automated design is in sharp contrast to current spacecraft design methods. Current design methods use engineering experience and judgement of a team of experts to identify candidate designs. This process requires significant expertise and experience and is long and expensive in terms of time and labor. The lack of a systematic approach to fully evaluate the whole design space might lead to a sub-optimal solution or worse an intractable solution. Such an approach can also miss innovative designs not thought of by the human design team. Evolutionary design techniques however benefit from the exponential rise in computational speed and can overcome human team limitations, including number of experts available, time available to perform a design task, and time available to meet together and work continuously. A machine-automated approach can augment a human design team by being able to review a larger candidate of designs than would otherwise be possible and thus find better designs.
Our approach to the spacecraft design problem is modelled after the knapsack problem (KP). Optimization of the knapsack problem is considered an NP (non-deterministic polynomial-time) hard problem. NP is the set of decision problems in computational complexity theory that are solvable by a non-deterministic Turing machine in polynomial time. NP hard are the problems that are at least as hard as the hardest problems in NP and they are not necessarily decision problems [
17]. In the knapsack problem, the goal is filling a knapsack with as many items so that the value of the items is maximized without exceeding the capacity of the knapsack. For a spacecraft, thanks to high launch costs, mass and volume are still a premium. In our spacecraft design problem, the focus is to effectively package components from a COTS inventory within a specified mass and volume constraint satisfying the standard CubeSat form factor. In our approach, the design of the spacecraft is specified by a gene that determines what components are to be packaged inside the spacecraft. This method of describing the design with genes keeps the process simple and produces results fast using desktop computers. In the following sections, we provide a brief description of the knapsack problem and the multidimensional multiple-choice knapsack problem which is the basis for our CubeSat design problem through trade space selection. Next, we present the optimization problem to design a 3U nadir-pointing, earth-observing CubeSat in low Earth orbit (LEO) followed by the models used for simulation. This is followed by the results and conclusions.
2. Knapsack Problem
The knapsack problem is analogous to the packing of items into a knapsack, where items are selected from an available set such that the total value of the selected items is maximized while simultaneously not exceeding the knapsack’s capacity [
18]. It is a highly popular problem in the research field of constrained and combinatorial optimization. The most common problem being solved is the 0-1 knapsack problem, which restricts the number of copies,
, of each kind of object,
, to zero or one. Given a set of
objects
, each with a weight,
, and a value,
, along with a maximum weight capacity,
, of the knapsack, the objective is to find a subset,
, in such a way that the weight sum over the objects in
does not exceed the knapsack capacity and yields a maximum value. The problem is mathematically expressed as Equation (1).
Here, represents the number of instances of item to include in the knapsack. Informally, the problem is to maximize the sum of the values of the items in the knapsack so that the sum of the weights is less than or equal to the knapsack’s capacity. The problem is a reduction of real-world resource allocation problems with constraints. The 0–1 KP is weakly NP-hard, which means it can be solved exactly by utilizing pseudo-polynomial time algorithms based on dynamic programming. There are also other algorithms that can solve it a very short amount of time such as greedy algorithms and core branch and bound algorithms. However, the difficulty of the problem increases as the number of variables and correlation between them is increased.
3. Knapsack Problem in System Design
The knapsack problem is also relevant for designing systems consisting of multiple subsystems such that its performance is maximized satisfying multiple resource constraints. Considering a system to be composed of
subsystems with each subsystem,
, consisting of
options, as shown in
Figure 1, each option of the subsystem has a particular value and it requires
resources. The objective is to pick exactly one option from each subsystem for the maximum total value of the collected items, subject to
resource constraints of the knapsack (system).
This is a variant of the classical 0-1 knapsack problem and is termed the multidimensional multiple-choice knapsack problem (MMKP) [
19]. In mathematical notation, let
and
be the value and resource vector of the object
, i.e.,
option of the
subsystem. In addition, assume that
is the resource bound to the knapsack (system). Now, the problem is mathematically expressed as Equation (2).
The number of possible designs,
, is the product of the number of options of all
subsystems;
, which grows exponentially as the number of subsystems,
, and number of options,
, increases as shown in
Figure 2a. One way to solve this problem is to use a brute force method by trying all
possible designs, which is highly inefficient. Other methods include greedy algorithms and dynamic programming [
18]. Several heuristic algorithms have also been proposed for MMKP. One such algorithm is based on the Lagrange multiplier method [
20]. Other methods include HEU [
21], the modified guided local search [
22], the modified reactive local search [
23], the convex hull-based method [
24], and the column generation method [
25]. Several exact algorithms have also been proposed such as the branch-and-bound method with linear programming [
21] and the EMKP algorithm [
26]. Although these algorithms are able to obtain fairly good solutions for small search spaces, they become ineffective due to the
curse of dimensionality when approaching practical scenarios, where the number of subsystems and the number of options for each subsystem grows significantly.
However, EA provides a way to solve the knapsack problem in linear time. Evolutionary algorithms gain their inspiration from the natural process of evolution and emulate evolution by applying recombination (crossover), mutation, and natural selection on a population [
27]. Highly fit individuals from the population are sampled, recombined, and resampled to form offspring of potentially higher fitness. As the evolutionary algorithm progresses through generations, fitter building blocks proliferate and unfit building blocks are culled from the population. In a way, by working with these particular schemata the complexity of the problem is reduced. Instead of building a high-performance individual by trying every conceivable combination, better and better individuals are constructed from the best partial solutions of past sampling.
Evolutionary algorithms (EAs) are a stochastic optimization method that provides an approach to learning by mimicking the metaphor of natural biological evolution or Darwinian evolution. By applying the principle of survival of the fittest, EAs progressively produce a better solution over generations by operating on a population of potential solutions. By mimicking the process of natural evolution, EAs perform evolutionary operations borrowed from natural genetics such as selection, crossover, and mutation. Theoretically, this method of evolving a population of individuals would progressively generate individuals that are better suited to their environment, just as in natural adaptation [
28,
29].
Figure 2b shows the structure of a simple evolutionary algorithm (EA). When successful, this approach will produce near-optimal solutions. These solutions when compared to solutions obtained through point design by a team of experts potentially presents major improvements. However, these solutions should be used as a starting point by the design team and would need further study to determine their implementation feasibility.
4. CubeSat Design
Current CubeSat design methods require a team of experts, who use their engineering expertise and decision to develop a spacecraft design from available COTS components. Such an approach may not result in an optimal design in terms of mass, volume, and power capabilities and can also miss innovative designs not thought of by the human design team. As such, the CubeSat design problem is modeled as a multidimensional multiple-choice knapsack problem. The goal is to effectively package subsystems from an inventory of COTS components within a specified mass and volume constraint such that the spacecraft can meet the defined mission goals.
As a test example, the objective is to design a 3U nadir-pointing, earth-observing CubeSat in LEO. The CubeSat is composed of many components that can be broken down into six major subsystems: (a) structure, (b) command and data handling (on-board computer), (c) communication (transceiver and antenna), (d) power (solar panels, battery pack, and power management board), (e) attitude determination and control unit (reaction wheels, magnetorquers, sun sensor, accelerometer, star tracker, magnetometer), and (f) payload (camera and lens) [
1]. At first, an inventory of all the components is created from major CubeSat vendors (GomSpace, ISIS, Nanoracks, EnduroSat, DHV, NanoAvionics, Maryland Aerospace, Blue Canyon, Clydespace, Pumpkin, Space Micro, Rincon Research) which is used as a database for the Knapsack Problem. To solve the MMKP, an EA is used where a candidate population is evolved for hundreds of generations using evolutionary operators like selection, crossover, and mutation. Each individual in the population is defined by a gene,
, which specifies what and how many components are to be packaged inside the spacecraft. This keeps the gene and design process simple and produces results fast using desktop computers. There are also other approaches to automated design, and these include the use of variable length generative coding schemes [
28] that generate a constructive program to design the gene. Other bio-inspired approaches model morphogenesis [
30].
An example gene structure is shown in
Table 1. It consists of two parts: the first part is the subsystem specifier {sID, cID, aID, tID, bID, pID, adID, spID, caID} that selects components from the available database, and the second part selects the number of batteries {
} and the number and location of solar panels: body panels {b1, b2, b3, b4}, side deployed panels {ds1, ds2, ds3, ds4}, and top deployed panels {dt1, dt2, dt3, dt4}.
Figure 3 shows an example of a phenotype of the CubeSat rendered in MATLAB VRML. A script for automated assembly of all the components inside the structure is written such that each component is modeled as a cube with its maximum dimensions, and a spacing of 5 mm between components is provided for the wiring space requirement.
The design objective is to maximize the total downloadable ground area coverage and minimize the total mass and total cost of the satellite. For a selected camera from the inventory, the spatial resolution,
, is calculated as Equation (3).
where
is the altitude of the satellite,
is the focal length of the lense, and
is the pixel pitch of the image sensor as shown in
Figure 4. Considering the sensor having
number of pixels, the total ground area covered per image is given by
and the total size of each image is
, where
is the size of each pixel. If
is the total downloadable data, then the total downloadable ground area coverage,
, is expressed as Equation (4).
The mass and cost of each component is added to determine the total mass,
, and total cost,
, of the CubeSat, and thus, the objective function for the MMKP is expressed as Equation (5).
To complete the knapsack problem, several constraints are added in terms of mass, volume, power and other properties such as the state of charge (SOC) of the battery subsystem, the pointing accuracy of the ADCS subsystem, the frequency and storage capacity of the on-board computer (OBC), and the bandwidth of the antenna and transceiver. The first constraint is such that the total mass of the CubeSat is less than the maximum allowable mass of a 3U CubeSat (4 kg). For volume constraint, three constraints are added such that after all the components are stacked inside the structure the dimensions along x, y, and z direction (
) are within the 3U CubeSat standard (10 × 10 × 34 cm) [
31]. For power, two constraints are added: the first is that the average solar power generated per orbit,
, should be greater than the total power requirement of the CubeSat,
, (which is determined from each component selected), and the second is that the state of charge of the battery subsystem should not drop below a minimum value (user-defined). Other constraints are also added: ADCS—pointing accuracy is within a user-defined value, OBC—clock frequency,
and storage capacity are within some user-defined values, transceiver and antenna—bandwidth of antenna selected matches that of transceiver. The complete MMKP CubeSat design problem is then mathematically expressed as Equation (6).
Now, to calculate the total downloadable data,
, average solar power generated,
, state of charge of the batteries,
, and pointing accuracy of the ADCS system, a simulation environment is created that models orbit dynamics, attitude dynamics, solar power generation, battery state of charge, and communication disciplines as shown in
Figure 5. Every individual of the population is an input to the COTS database that produces the CubeSat design. The design is then an input to the simulation environment, and then the outputs of the simulation environment along with mass, volume, power, and other properties of the design are used to determine the fitness (combination of objective function and the constraints). Details of the modeled disciplines are provided below.
4.1. Orbit Dynamics
The orbit-dynamics discipline computes the Earth-to-body position vector of the satellite in the Earth-centered-inertial (ECI) frame. In the orbit equation, the
,
, and
coefficients capture the fact that the Earth is not a perfectly spherical and homogeneous mass and is expressed by Equation (7).
where
is the position vector of the satellite in the ECI frame,
is Earth’s gravitational parameter,
is Earth’s radius and
,
, and
are orbit perturbation coefficients. The
,
, and
terms must be considered because their effect is to rotate the orbit plane, which affects the power generation and communication discipline.
4.2. Attitude Dynamics
Since the CubeSat is an earth observing satellite, it must always maintain a nadir-pointing orientation. The attitude dynamics is simulated based on the reaction wheels/magnetorquers available in the selected ADCS system from the COTS inventory according to Equation (8).
where
is the moment of inertia matrix of the satellite,
is the angular velocity vector of the satellite,
is the moment of inertia matrix of the reaction wheel system,
is the angular velocity vector of the reaction wheel system,
is the torque applied by the reaction wheels system,
is the torque applied by the magnetorquers, and
is the external disturbance torques, which is a sum of the gravity gradient, aerodynamic, and solar radiation pressure disturbances. The simulation of the attitude dynamics based on the selected ADCS system provides the pointing accuracy of the satellite, which is used for the optimization problem.
4.3. Solar Power Generation
For solar power, first the sun-to-satellite line-of-sight,
, is calculated which is essentially a multiplier for the solar panel exposed area to calculate the solar power generated. The
is 0 if the satellite is in the eclipse and 1 otherwise, however, to implement the umbra and penumbra effects, the discontinuous function is smoothed by using a cubic function [
32], as shown in Equation (9).
where
and
are given by Equation (10).
The vectors
and
are shown in
Figure 6. The value of
represents how far this smoothing effect extends into the Earth’s shadow, typically
. Next, with the position of the sun known with respect to the satellite, the exposed solar panel area is required to calculate the total solar power generated. The exposed solar panel area is done in two steps. First for each panel having an area
and unit normal vector
and the body-to-sun unit vector
, the projected area of each panel on a plane perpendicular to
is given by Equation (11) and shown in
Figure 7a.
Next, if the satellite consists of deployable solar panels, the effects of shadows projected by one panel onto another must be considered. For two panels
with areas
and unit normal vectors
, and body-to-sun unit vector
, the shadow cast by
on
is calculated by two successive projections and an intersection operator as shown in
Figure 7b. The first projection is of
on a plane perpendicular to
given by
, and the second projection is of
on a plane perpendicular to
given by
. Finally, the projected shadow cast by
on
on a plane perpendicular to
is the intersection of areas
and
, as shown by Equation (12).
Thus, the total exposed solar panel area is expressed as Equation (13).
where
is the total number of solar panels, and
are the neighboring panels for panel
. Finally, the solar power generated is calculated as Equation (14).
where
is the solar constant, and
is the efficiency of the photovoltaic solar cells on each solar panel.
Figure 8 shows the simulation results of the solar panel exposed area model for a CubeSat with two panels: one body panel and one deployed panel, as shown in
Figure 8a. The simulation is performed in MATLAB, and the rendering of the CubeSat model shown in
Figure 8a is done in MATLAB VRML. The body-to-sun unit vector
is represented in a spherical coordinate system with respect to the body frame, and as such, the results are presented against azimuth and elevation angles of the sun.
Figure 8b shows the exposed area of panel 1 ignoring shadows cast by panel 2, while
Figure 8c shows the projected shadow cast by panel 2 on panel 1, and finally,
Figure 8d shows the total exposed area of panel 1. Our method is able to calculate the total exposed solar panel area for the above CubeSat with two panels, which is used for more complicated designs evolved through the optimization process.
4.4. Battery State of Charge
The state of charge of the battery system is determined by the Ordinary Differential Equation (ODE) (Equation (15)) [
32].
where
is the nominal discharge capacity of the battery system with
being the number of batteries and
being the discharge capacity of each battery. The dependence of the battery voltage is on the SOC and also on the temperature, which is exponential, as shown by Equation (16).
where
is the reference temperature, and
is the temperature decay coefficient. At any given time, the battery power is the sum of the loads, given by Equation (17).
where
is the power generated by the solar panels,
is the power consumed by the ADCS system,
is the communication power, and
is the power consumed during the nominal operation of the CubeSat and the payload (camera).
4.5. Communication
For communication, first the ground station-to-satellite line-of-sight,
, is computed based on the dot product between the Earth-to-ground station unit vector and the ground station-to-body vector in the inertial frame, as shown in
Figure 9. The discontinuous function is smoothed using a cubic function. Next, the resulting data download rate,
, is calculated using Equation (18) [
32,
33].
where
is the speed of light,
is the Boltzmann constant,
is the transmission frequency,
is the system noise temperature,
is the receiver (ground station) gain,
is the maximum acceptable signal to noise ratio,
is the communication efficiency,
is the line loss factor,
is the communication power,
is the transmitter (satellite) gain, and
is the distance from the satellite to the ground station.
4.6. Optimization
Now to solve the described CubeSat design problem implemented as a multidimensional multiple-choice knapsack problem, an evolutionary algorithm is used. The constraints described in Equation (6) divide the search space into two regions—feasible and infeasible regions. The constraints are handled using the penalty function approach. For each solution
,
being the index of the gene, the constraint violation for the inequality constraints
are calculated as in Equation (19).
Thereafter, all constraint violations are added together to get the overall constraint violation as in Equation (20). There are 10 constraints in the optimization problem and constraint 10 is decomposed into 2 inequality constraints.
This constraint violation is then multiplied with a penalty parameter,
, and then the product is added to each of the objective function as in Equation (21).
The fitness function,
, takes into account the constraint violations. For a feasible solution the corresponding
term is zero and
becomes equal to the original objective function
. However, for an infeasible solution,
, thereby adding a penalty corresponding to the total constraint violation. Once the fitness function is derived, the steps for an EA are used to optimize the problem. Initially, a random parent population,
, is created of size
. Each individual of the parent population comprises of the design variables represented by the gene,
, chosen by a discrete uniform distribution
, where
is the j-th design variable and
and
are the upper and lower bounds of the design variables. For each individual, the value of the fitness function
is evaluated using the simulation environment as shown in
Figure 6. Next, the usual tournament selection, crossover, and mutation operators are used to create an offspring population,
, of size
[
34]. For crossover, a blend crossover (BLX-α) operator is used [
35]. During mutation, we use a non-uniform mutation operator. The mutation operator is constructed such that the probability of creating a solution closer to the parent is more than the probability of creating one away from it [
36]. By comparing the current population with previously found best solutions, elitism is introduced into the algorithm. For the
generation, first a combined population
is formed. The population
is of size
. Since all previous and current population members are included in
, elitism is automatically ensured. Next, the fitness of each individual in
is calculated according to the fitness function
. Based on the value of the fitness function, the population is sorted in ascending order and the best
individuals are selected for the next generation. The new population,
, of size
is again used for selection, crossover, and mutation to create a new population,
, of size
. The process is repeated until the desired results are obtained, or the max number of generations are achieved.
5. Results
For the test case the orbit of the satellite is simulated for an altitude of 400 km, inclination of 90°, eccentricity of 0 (circular orbit), RAAN of 90°, argument of perigee of 0°, and mean anomaly of 0°. The values of the constants used in the simulation are shown in
Table 2. The optimization problem is simulated for 10 orbits in MATLAB, and the results are presented in
Figure 10 and
Figure 11.
Figure 10a shows the evolution of the fitness function
over 250 generations with 100 populations in each generation.
Figure 10b shows the evolution of the total downloadable ground area coverage,
,
Figure 10c shows that of the total mass,
, and
Figure 10d shows that of the total cost,
, of the CubeSat over 250 generations. The shaded region of each of the plots show the 95% confidence interval over 20 runs. The results of the plots show that the fitness function decreased with generations converging to a constant fitness, while the total downloadable ground area coverage increased, and the total mass and total cost decreased with generations as implemented in the optimization problem.
Figure 11 shows the evolution of the constraints
to
. The constraint
is already presented in
Figure 10c, which shows that the total mass is below the 4 kg allotted mass of a 3U CubeSat.
Figure 11a,b shows the evolution of the dimension of the assembled CubeSat components along z, x, and y directions, showing it is within the 340 mm and 105 mm limit (shown by the dashed lines).
Figure 11c shows the evolution of the average solar power generated per orbit,
, and the total power requirement of the CubeSat,
, satisfying the requirement
.
Figure 11d shows the evolution of the state of charge, SOC, of the battery system with the dashed line showing the minimum requirement of 0.2.
Figure 11e shows the evolution of the pointing accuracy of the ADCS system, which satisfies the constraint
.
Figure 11f,g shows the evolution of the clock frequency and storage capacity of the on-board computer, which are greater than the specified values of 730 MHz and 1.25 Gb, respectively (shown by the dashed lines).
Figure 11h shows the evolution of the antenna frequency (
) along with the upper (
) and lower (
) limit of the frequency of the transceiver, showing the antenna frequency is within the limit of the bandwidth of the transceiver. It is clearly evident that all the constraints have been satisfied when the optimization process was terminated, and as such, the final solutions found are in the feasible region of the design space.
Figure 12 shows the rendering of the evolution of the best design over generations 1, 12, 15, 20, 30, 80, 100, and 150. It can be seen that the number of solar panels and their positions changed over generations, and at the same time, the internal components vertically stacked inside the structure also changed.
Figure 13 and
Figure 14 show the variation of the design variables (represented by the gene in
Figure 3) of the best individual over generations. The color of the design variables represents the component selected from the inventory in
Figure 13 and the number of batteries and solar panels in
Figure 14. It can be seen that the combination of the subsystems and the number of batteries and solar panels changed until generation 150, and then an optimal combination was found.
Finally, we present a comparative analysis between three stochastic algorithms, evolutionary algorithm (EA), particle swarm optimization (PSO) [
37], and simulated annealing (SA) [
38], in terms of performance on the CubeSat design MMKP problem. The theory behind the application of the EA to our design problem has already been discussed in the preceding sections. We now provide a brief explanation on how PSO and SA are implemented.
5.1. Particle Swarm Optimization
Particle swarm optimization is a stochastic optimization method that is inspired from the social behavior of fish schooling or birds flocking. Just like EA, PSO is also initialized with a population (swarm) of individuals (particles). However, unlike EA, PSO does not use evolutionary operations to search for optima; instead, each particle is guided through the search space in every iteration by updating their positions and velocities using better positions found by other particles. In our CubeSat design problem, the fitness is defined by Equation (21) as
, where
is the position of each particle, and the velocity of each particle is denoted by
. The algorithm is initialized by generating
particles with their positions and velocities initialized randomly. Next, at every iteration,
, the velocity and position of each particle are updated by Equations (22) and (23), respectively.
where
,
, and
are the tuning parameters,
and
are random numbers generated between 0 and 1,
is the best remembered individual particle position, and
is the best remembered swarm position. This process of updating the position and velocity of each particle continues until the desired fitness is achieved or maximum number of iterations reached.
5.2. Simulated Annealing
Simulated annealing is another stochastic optimization method that is inspired from annealing in metallurgy. The algorithm starts with a random trial solution,
, whose fitness is determined by
, defined by Equation (21), and a large initial temperature,
. At each iteration,
, the algorithm randomly generates a new solution,
, and its corresponding fitness,
, is calculated. If the new solution is better than the current solution, it is automatically selected as the next solution. However, if the new solution is worse, the algorithm accepts it as the next solution based on a probability of acceptance function defined by Equation (24).
where
and
is the current temperature, which is systematically lowered by the algorithm at each iteration according to Equation (25).
where
is the annealing parameter. Since the values of Δ and
are positive, larger Δ and a lower temperature leads to a smaller acceptance probability. The algorithm stores the best solution found so far and proceeds until the desired fitness is achieved or maximum number of iterations reached.
Figure 15a shows the evolution of the mean and standard deviation of fitness function over generations for all the 3 algorithms simulated and averaged over 20 times. The constant parameters used for PSO and SA are presented in
Table 2. It can be clearly seen that the optimal solution found using simulated annealing lies in the infeasible region as the fitness function is positive. Both the evolutionary algorithm and particle swarm optimization were able to find a feasible solution, but the performance of the EA is better than PSO, as shown in
Figure 15b, making it the better choice for this problem. Our observations show PSO stagnates prematurely, while EAs are more effective in continuing to innovate over generations.