The optimization of crop rotations has a long history in agricultural operations research. Early work employed Linear Programming (LP) to allocate crops to fields and seasons based on profit maximization subject to resource constraints [23], [24]. Detlefsen and Jensen [13] extended LP-based approaches to incorporate crop sequence restrictions, demonstrating feasibility for Danish farming systems but acknowledging limitations in capturing nonlinear soil dynamics and stochastic weather effects. MILP formulations have been applied to crop rotation problems with combinatorial crop assignment decisions [14], [25], though computational tractability becomes problematic as the number of parcels, years, and candidate crops increases.

Rule-based approaches, drawing on agronomic expert knowledge, have been widely used in European farming systems. ROTOR [11] implemented a knowledge-based system for generating crop rotations satisfying ecological and agronomic constraints in organic farming. While effective for constraint satisfaction, such systems do not optimize across continuous objective spaces and struggle to incorporate climate uncertainty. The CropRota model [26] combined LP with agronomic rules for Austrian agriculture but was limited to a single economic objective.

Evolutionary computation has brought new capabilities to crop planning. Annetts and Audsley [15] applied GAs to whole-farm planning, including crop selection, in UK conditions. Dury et al. [10] used a multi-objective GA for crop allocation in French farming systems, demonstrating the feasibility of Pareto-based exploration but without climate integration. More recently, Akplogan et al. [27] employed constraint programming combined with metaheuristics for vegetable crop rotation under agronomic constraints. Kosmidou [28] explored GA-based crop pattern optimization for water-scarce regions. Reinforcement learning approaches have also emerged, with Gautron et al. [29] applying deep reinforcement learning to crop management decisions, though not specifically to multi-year rotation sequencing under climate uncertainty.

MOEAs have been applied across a range of agricultural planning problems. NSGA-II [16] has been particularly prevalent, employed in irrigation scheduling [30], [31], land-use allocation [19], and water resource management [32]. MOEA/D [17] has seen application in resource allocation problems with many objectives, though its use in agriculture remains less common than NSGA-II. SPEA2 [33] has been employed in agricultural water management and crop planning contexts.

MOPSO, introduced by Coello Coello et al. [20] using an external archive and adaptive grid mechanism, has been applied to engineering design optimization and resource scheduling. In agricultural contexts, Nguyen et al. [34] applied MOPSO to precision irrigation scheduling, while Sarker and Ray [35] explored swarm-based methods for crop planning in Australian conditions. Memmah et al. [18] used multi-objective metaheuristics for landscape-level crop allocation in France, integrating biodiversity and economic objectives. However, none of these MOPSO applications have addressed multi-year crop rotation planning with integrated climate projections and soil dynamics.

The adaptation of swarm intelligence to discrete–combinatorial search spaces, as required for crop assignment decisions, has been explored through binary PSO [36], discrete PSO variants [37], and Hybrid approaches combining PSO velocity mechanics with permutation-based operators [38]. These techniques provide a foundation for the discrete–continuous hybrid encoding required in crop rotation optimization, where discrete crop assignments must be jointly optimized with continuous agronomic parameters.

The integration of climate projections into agricultural decision-making has accelerated with the availability of CMIP5 and CMIP6 ensemble outputs [39], [40]. Scenario-based optimization under climate uncertainty has been explored using robust optimization frameworks [41], stochastic programming [42], and scenario sampling approaches [43]. In agricultural planning, Thornton et al. [44] and Challinor et al. [45] have assessed climate impacts on crop production using ensemble crop model simulations driven by multiple GCM outputs, establishing the importance of multi-model, multi-scenario assessment.

Process-based crop models serve as essential tools for translating climate inputs into yield responses. Decision Support System for Agrotechnology Transfer (DSSAT) [46] provides detailed simulation of crop growth, soil water, and nutrient dynamics for over 40 crops. Agricultural Production Systems Simulator (APSIM) [47] offers modular simulation capabilities widely applied in tropical and dryland agriculture. AquaCrop [48], developed by FAO, emphasizes water-limited yield prediction with a relatively parsimonious parameter set, making it suitable for data-scarce environments. Sultan et al. [8] applied CMIP6 projections to assess climate change impacts on West African cereal yields, finding that CMIP6 models project less negative impacts than CMIP5, highlighting the importance of using current-generation climate models in agricultural planning. Despite these advances, the direct integration of climate scenario ensembles into multi-objective optimization algorithms for crop rotation planning has not been systematically explored. Most existing studies either use climate projections to assess pre-defined rotation scenarios ex post or optimize under a single deterministic climate assumption, neglecting the distributional uncertainty inherent in ensemble projections.

Soil nutrient dynamics are central to evaluating rotation quality but are frequently oversimplified in optimization frameworks. The CENTURY model [49] simulates long-term soil organic matter dynamics and nutrient cycling across multiple soil pools. The RothC model [50] focuses on turnover of soil organic carbon. These process-based soil models capture the accumulation and depletion dynamics relevant to rotation evaluation but are computationally expensive when embedded within iterative optimization loops requiring thousands of fitness evaluations.

Simplified soil nutrient accounting models have been developed for integration with optimization frameworks. Dogliotti et al. [12] incorporated nitrogen balance calculations into a rotation design tool for Uruguayan farming systems. Reckling et al. [4] used simplified nitrogen and phosphorus budgets to evaluate the environmental performance of diversified rotations across European sites. These approaches sacrifice mechanistic detail for computational tractability, a trade-off that is particularly relevant in metaheuristic optimization where population-based search requires large numbers of candidate evaluations.

The explicit incorporation of dynamic soil state variables soil organic matter, available nitrogen, available phosphorus, and soil pH that evolve across rotation years as a function of crop-specific uptake, residue inputs, and management practices has not been embedded within a multi-objective swarm optimization framework for crop rotation planning. This integration is essential for evaluating the long-term sustainability of candidate rotations, as rotations that appear economically optimal in the short term may drive cumulative soil degradation over multi-year horizons.

We formulate the climate-adaptive crop rotation planning problem as a constrained multi-objective optimization problem. Consider a farm comprising P parcels over a planning horizon of T years. Let S = {s1, s2, …, sM} denote the set of M candidate crops (including fallow as a pseudo-crop). The primary decision variable is the crop assignment matrix C = [C ] ∈ ST×P, where C ∈ S specifies the crop assigned to parcel p in

tp tp

year t. Additionally, continuous decision variables u = (u irr, u fert, u seed) represent irrigation volume,

tp tp tp tp

fertilizer application rate, and seeding density for each crop–parcel–year combination, bounded within agronomically feasible ranges.

The problem seeks to simultaneously optimize three objective functions:

Objective 1. Maximize expected economic return (f1):

f (C, u) = (1/K) ∑ K ∑ T ∑ P [Y(C , θ k, σ ) × π(C ) − Γ(C , u )] × (1+r)−t,

1 k=1 t=1 p=1 tp tp tp tp tp tp (1)

where Y(·) is the crop yield function dependent on crop type C , climate variables θ k under scenario k, and

tp tp

soil state σtp; π(·) is the market price per unit yield; Γ(·) is the total production cost including inputs, labor, and mechanization; r is the discount rate; and K is the number of climate scenarios. The yield function is detailed in Section 3.3.

Objective 2. Minimize environmental degradation (f2):

f2(C, u) = α1 ΔNcum + α2 Wcum + α3 ECO2, (2)

where ΔNcum = ∑t,p max(0, Nuptake(Ctp) − Ninput(Ctp, utp) − Nfix(Ctp)) is the cumulative soil nitrogen depletion; Wcum is the cumulative water consumption index normalized by regional water availability; ECO2 is the cumulative carbon emission index from tillage, fertilizer application, and fuel use; and α1, α2, α3 are normalization weights ensuring commensurate scaling of the environmental sub-indices.

Objective 3. Maximize climate resilience (f3):

f (C, u) = min { ∑ T ∑ P Y(C , θ k, σ ) / Y (C ) },

3 k=1,…,K t=1 p=1 tp tp tp ref tp (3)

where Yref(·) denotes the reference yield under historical average climate conditions. This minimax formulation ensures that the optimizer explicitly seeks rotation plans whose performance degrades least under the most adverse climate scenario, providing a measure of distributionally robust planning.

The optimization is subject to the following constraint classes:

I. Crop succession constraints: no crop from the same botanical family may be planted on the same parcel in consecutive years: fam(Ctp) ≠ fam(C(t+1)p) for all t, p.

II. Water budget constraints: total irrigation demand across all parcels must not exceed seasonal regional water availability: ∑ u irr ≤ W (t) for all t.

p tp avail

III. Labor constraints: total labor demand must not exceed household or seasonal labor supply: ∑p L(Ctp) ≤ Lmax(t) for all t.

IV. Minimum dietary diversity: at least Dmin distinct food groups must be represented in each year's crop portfolio: |{food_group(Ctp) : p = 1,…,P}| ≥ Dmin for all t.

V. Market absorption limits: total production of any single crop across all parcels must not exceed the local market absorption threshold: ∑p Y(Ctp,·) ≤ Qmax(Ctp) for all t.

VI. Fallow minimum: at least one parcel per rotation cycle must undergo fallow for soil recovery.

MOSWARM-CRP extends the canonical MOPSO framework with domain-specific innovations tailored to the crop rotation search space. A swarm of N particles collectively searches the space of feasible rotation plans. Each particle i encodes a candidate solution (Ci, ui) and maintains a velocity vector governing its trajectory through the search space.

Particle encoding: the hybrid discrete–continuous encoding represents each particle as a concatenation of two segments: 1) a discrete segment encoding the crop assignment matrix C ∈ ST×P, where each gene takes an integer value from {1, 2, …, M} corresponding to a candidate crop, and 2) a continuous segment encoding the agronomic input parameters u ∈ ℝT×P×3 (irrigation, fertilizer, seeding density), bounded within crop-specific feasible ranges. The total dimensionality of each particle is T × P × (1 + 3) = 4TP.

Velocity update: For the continuous segment, the standard PSO velocity update is applied with modifications:

v cont(t+1) = w(t) · v cont(t) + c r [pbest cont − x cont(t)] + c r [gbest cont − x cont(t)],

id id 1 1 id id 2 2 id id (4)

where w(t) is the adaptive inertia weight (Section 3.5), c1 and c2 are cognitive and social acceleration coefficients, r1 and r2 are uniform random numbers in [0,1], pbestid is the personal best position of particle i in dimension d, and gbestid is a global guide selected from the Pareto archive using crowding distance-based tournament selection.

For the discrete segment (crop assignments), a probabilistic swap operator is employed. At each dimension d corresponding to a crop assignment Ctp, the crop value is retained, swapped to the personal best value, or swapped to the archive guide value according to probabilities derived from the velocity magnitude:

P (d) = sigmoid(v disc) = 1 / (1 + exp(−v disc)),

swap id id (5)

where the discrete velocity v disc is updated analogously to the continuous case but operates in the

id

probabilistic swap space, following the binary PSO formulation of Kennedy and Eberhart [36] extended to the multi-valued discrete case.

Pareto archive management: an external archive A of maximum size NA stores non-dominated solutions discovered during the search. When the archive exceeds capacity, solutions in the most crowded regions of the objective space are pruned using an adaptive grid mechanism [20]. The objective space is partitioned into hypercubes, and solutions residing in hypercubes with the highest occupancy are preferentially removed. Crowding distance [16] is computed within each hypercube to preserve boundary and extreme solutions. Archive guide selection: when selecting a global guide gbest from the archive for each particle, a two-stage tournament is employed: 1) a hypercube is selected with probability inversely proportional to its occupancy (favoring guides from sparse regions), and 2) within the selected hypercube, the solution with the largest crowding distance is chosen. This mechanism promotes exploration of under-represented Pareto front regions.

The main algorithmic loop of MOSWARM-CRP is presented as Algorithm 1:

Algorithm 1. MOSWARM-CRP main loop.

Input: Swarm size N, archive size N , max iterations T , climate scenarios {θ1,…,θK}, crop set S, soil parameters,

A max

constraints Output: Pareto-optimal archive A of rotation plans 1: Initialize swarm: for i = 1 to N do 2: Generate random feasible rotation plan (Ci, ui) 3: Apply constraint repair operator (Section 3.6) 4: Initialize velocity vi ← 0 5: Set pbesti ← (Ci, ui) 6: end for 7: Evaluate all particles: for each i, compute (f , f , f ) across all K climate scenarios

1 2 3

8: Initialize archive A with non-dominated solutions 9: Construct adaptive grid over A 10: for t = 1 to Tmax do 11: Compute adaptive inertia weight w(t) (Algorithm 2) 12: for i = 1 to N do 13: Select guide gbesti from A via crowding tournament 14: Update continuous velocity (Eq. 4) 15: Update discrete velocity and apply swap (Eq. 5) 16: Update position (Ci, ui) 17: Apply constraint repair operator (Section 3.6) 18: Evaluate (f , f , f ) across K scenarios 19:

1 2 3

Update pbesti using Pareto dominance 20: end for 21: Update archive A with new non-dominated solutions 22: If |A| > NA: prune via adaptive grid crowding 23: Update adaptive grid boundaries 24: Record convergence metrics (HV, archive diversity) 25: end for 26: return A.

A distinctive feature of MOSWARM-CRP is the explicit evaluation of each candidate rotation plan across an ensemble of K = 10 climate scenarios derived from CMIP6 [39]. Five Global Climate Models (GCMs), CanESM5, MIROC6, MPI-ESM1-2-HR, GFDL-ESM4, and UKESM1-0-LL, are employed under two SSP pathways (SSP2-4.5 representing moderate mitigation and SSP5-8.5 representing high-emission futures), yielding 5 × 2 = 10 distinct climate trajectories for the planning period 2025–2050. Raw GCM outputs are statistically downscaled to the study regions using quantile mapping [51] calibrated against ERA5 reanalysis data [52] and local meteorological station records.

For each climate scenario k and each year–parcel combination (t, p), the relevant climate variables are extracted: daily maximum and minimum temperature (Tmax, Tmin), daily Precipitation (P), and daily solar Radiation (Rs). These are aggregated into growing-season summary variables used by the yield function. Crop yield response model: we employ a simplified, AquaCrop-inspired yield response function [48], [53] that relates crop yield to water stress and temperature stress through Growing Degree Days (GDD) accumulation:

Y(s, θk, σ) = Y (s) × K (s, θk) × K (s, θk) × K (s, σ) × K (s, θk),

pot w T soil GDD (6)

where Ypot(s) is the potential yield of crop s under non-stressed conditions (derived from FAO and regional statistics); Kw ∈ [0,1] is the water stress coefficient computed from the ratio of actual to potential evapotranspiration; KT ∈ [0,1] is the temperature stress coefficient capturing heat and cold stress effects; Ksoil ∈ [0,1] is the soil fertility coefficient dependent on available nitrogen and phosphorus relative to crop requirements; and KGDD ∈ [0,1] penalizes incomplete thermal accumulation when GDD falls below the crop's maturity requirement.

The water stress coefficient is calculated as

Kw = 1 − ky(s) × (1 − ETa / ETc),

where ky(s) is the FAO yield response factor for crop [54], ETa is actual evapotranspiration estimated from precipitation and irrigation inputs, and ETc is crop potential evapotranspiration computed using the Penman-Monteith equation with crop coefficients [55]. The temperature stress coefficient is modeled as a trapezoidal function of mean growing-season temperature relative to crop-specific cardinal temperatures (Tbase, Topt,low, Topt,high, Tmax).

An embedded soil nutrient dynamics model tracks the evolution of four soil state variables across rotation years for each parcel: soil organic matter (SOM, in g/kg), available nitrogen (Navail, in kg/ha), available phosphorus (Pavail, in kg/ha), and soil pH. The model is initialized with parcel-specific baseline soil conditions and updated annually based on crop-specific nutrient dynamics.

Nitrogen balance: The annual nitrogen balance for parcel p in year t is:

Navail(t+1, p) = Navail(t, p) + Nfert(t, p) + Nfix(Ctp) + Nmin(SOMtp) + Nres(Ctp) − Nuptake(Ctp, Ytp) − Nleach(t, p) − Nvol(t, p),

where Nfert is applied fertilizer nitrogen; Nfix is biological nitrogen fixation (positive only for leguminous crops: groundnut, cowpea, chickpea, lentil, faba bean); Nmin is mineralization from soil organic matter (modeled as a fraction of SOM with temperature-dependent rate); Nres is nitrogen returned through crop residue incorporation; Nuptake is crop nitrogen uptake proportional to yield; Nleach is leaching losses estimated from precipitation surplus; and Nvol is volatilization losses.

Soil organic matter dynamics: SOM evolves through residue inputs and decomposition:

SOM(t+1, p) = SOM(t, p) × (1 − kdec) + h × Rcrop(Ctp, Ytp),

where k is the annual decomposition rate (climate- and texture-dependent, typically 0.02–0.05 yr−1), h is the

dec

humification coefficient (fraction of residue carbon converted to stable SOM), and Rcrop is the crop residue biomass returned to the soil. Leguminous crops and high-biomass cereals (sorghum, millet) contribute more residue than low-biomass crops, creating positive SOM trajectories under diversified rotations. Fallow periods allow enhanced mineralization and SOM equilibration.

Phosphorus and pH dynamics: available phosphorus is modeled with a simplified balance analogous to nitrogen, accounting for fertilizer inputs, crop uptake (proportional to yield and crop-specific P demand), and slow release from soil mineral P reserves. Soil pH evolves incrementally based on crop-specific rhizosphere effects, fertilizer type (ammonium-based fertilizers cause acidification), and lime application when included in the management variable set.

The soil model is computationally lightweight (requiring only arithmetic operations per parcel per year) yet captures the essential inter-annual feedback dynamics that determine whether a rotation sequence is building or depleting soil capital over the planning horizon.

The inertia weight w(t) in the velocity update (see Eq. (4)) governs the balance between exploration (high w) and exploitation (low w). Standard MOPSO implementations use either a linearly decreasing schedule [56] or a fixed value. MOSWARM-CRP introduces an adaptive inertia weight that responds to three feedback signals specific to the agricultural optimization domain:

w(t) = wbase(t) + β1 ΔHV(t) + β2 Darchive(t) + β3 φseason(t), (10)

where wbase(t) = wmax − (wmax − wmin) × (t / Tmax) is the linearly decreasing baseline (with wmax = 0.9, wmin = 0.4); ΔHV(t) is the normalized change in hypervolume indicator over the last τ = 20 iterations (positive when the Pareto front is improving, triggering exploitation; negative when stagnating, triggering exploration); Darchive(t) is an archive diversity metric based on the coefficient of variation of crowding distances in the archive (high diversity encourages exploitation, low diversity encourages exploration); φseason(t) is a seasonal agricultural calendar term that modulates exploration intensity to focus search on rotation configurations that respect planting and harvesting windows for the study region's agricultural calendar; and β1, β2, β3 are scaling parameters calibrated via preliminary experiments (set to 0.05, 0.03, and 0.02, respectively).

Algorithm 2. Adaptive inertia weight computation.

Input: Iteration t, archive A, HV history, calendar constraints Output: Inertia weight w(t) 1: Compute baseline: wbase ← wmax − (wmax − wmin) × (t/Tmax) 2: Compute HV improvement rate: 3: ΔHV ← (HV(t) − HV(t−τ)) / max(HV(t−τ), ε) 4: If ΔHV > 0: adj1 ← −β1 × ΔHV // exploiting 5: Else: adj1 ← β1 × |ΔHV| // exploring 6: Compute archive diversity: 7: CD ← crowding distances of all solutions in A 8: Darchive ← std(CD) / mean(CD) 9: If Darchive > δthresh: adj2 ← −β2 // diverse, exploit 10: Else: adj2 ← β2 // clustered, explore 11: Compute seasonal adjustment: 12: φseason ← penalty for off-calendar crop assignments 13: adj3 ← β3 × φseason 14: Combine: w(t) ← clip(wbase + adj1 + adj2 + adj3, wmin, wmax) 15: return w(t).

Constraint handling in MOSWARM-CRP employs a hybrid approach combining a repair operator with an adaptive penalty function, reflecting the heterogeneous nature of the constraints.

Repair operator: after each position update, the repair operator enforces hard constraints that define agronomic feasibility:

I. Crop succession repair: for each parcel, if consecutive years violate the same-family constraint, the later year's crop is replaced by the nearest feasible crop in the candidate set (selected by cycling through S in a deterministic order that preserves solution diversity).

II. Bound repair: continuous variables (irrigation, fertilizer, seeding density) exceeding their crop-specific bounds are clipped to the nearest feasible value.

III. Fallow enforcement: if no parcel in the rotation includes a fallow period within any consecutive Tfallow-year window, the lowest-profit parcel–year combination is replaced with fallow.

IV. Penalty function: soft constraints (water budget, labor, dietary diversity, market absorption) are handled through an adaptive penalty added to objective f2 (environmental degradation):

f penalized = f + λ(t) × ∑ max(0, g (C, u))2,

2 2 j j (11)

where gj(·) represents the j-th constraint violation magnitude and λ(t) is an iteration-dependent penalty coefficient that increases over the search to ensure convergence toward feasibility: λ(t) = λ × (1 + t / T )2.

0 max

Two case study regions were selected to represent contrasting agro-climatic zones in Africa, both characterized by high climate vulnerability and critical agricultural importance:

Ségou Region, Mali (Sahel, Semi-Arid Tropical): located in central Mali (13.4°N, 6.3°W), the Ségou Region lies within the Sudano-Sahelian zone, characterized by a unimodal rainy season (June–October) with mean annual precipitation of approximately 600 mm (range: 400–800 mm over 1980–2020). Agriculture is predominantly rainfed, with millet (pennisetum glaucum) and sorghum (sorghum bicolor) as staple cereals, complemented by groundnut (arachis hypogaea), cowpea (vigna unguiculata), and cotton (gossypium hirsutum) as cash or rotation crops. Soils are predominantly sandy (arenosols and lixisols) with low inherent fertility (SOM < 1%, total N < 0.05%). Smallholder farms of 3–10 ha dominate, with limited access to mineral fertilizers, mechanization, and irrigation. The region is a critical food security zone for Mali, supplying a significant share of national cereal production, and is highly exposed to rainfall variability and progressive warming trends.

Meknès-Fès Region, Morocco (Semi-Arid Mediterranean): Situated in north-central Morocco (33.9°N, 5.5°W), the Meknès-Fès Region experiences a Mediterranean climate with cool, wet winters and hot, dry summers, receiving approximately 450 mm mean annual precipitation concentrated between October and April. Agriculture is more diversified and intensified than in the Sahel, featuring durum wheat (Triticum durum), soft wheat (Triticum aestivum), and barley (Hordeum vulgare) as principal cereals, with food legumes (chickpea, lentil, faba bean), sunflower, and sugar beet as rotation crops. Supplemental irrigation is available for approximately 30% of cropland from reservoir systems. Soils are predominantly calcareous Vertisols and Cambisols with moderate fertility. Farm sizes range from 5–50 ha, and semi-commercial farming is common. The region faces progressive aridification, declining groundwater tables, and increasing inter-annual precipitation variability under climate change projections.

The candidate crop sets for each region, along with their key agronomic parameters, are presented in Tables 1 and 2. Parameters were sourced from FAO crop databases [7], regional agricultural statistics from the Institut d'Économie Rurale (IER) of Mali and the Institut National de la Recherche Agronomique (INRA) of Morocco, and published agronomic literature.

Historical data: daily meteorological records for the period 1980–2020 were obtained from national meteorological services (Direction Nationale de la Météorologie, Mali; Direction Générale de la Météorologie, Morocco) for stations nearest the study sites, supplemented by ERA5 reanalysis data [52] for gap-filling and spatial coverage extension. Variables include daily maximum/minimum temperature, precipitation, solar radiation, wind speed, and relative humidity.

Future projections: CMIP6 model outputs for the period 2025–2050 were downloaded from the Earth System Grid Federation (ESGF) data nodes. Five GCMs were selected based on demonstrated skill in reproducing African climate variability [57], [58]:

Statistical downscaling was performed using quantile delta mapping [59] to adjust daily GCM outputs to the local observational climate, preserving projected change signals while correcting distributional biases. Each GCM was run under SSP2-4.5 and SSP5-8.5, yielding K = 10 climate scenarios that span a plausible range of future climatic conditions for each study region.

MOSWARM-CRP was benchmarked against six alternative methods representing diverse optimization paradigms:

I. NSGA-II [16]: the canonical multi-objective GA with non-dominated sorting and crowding distance.

Population size 200, SBX crossover (pc = 0.9, ηc = 20), polynomial mutation (pm = 1/n, ηm = 20).

II. MOEA/D [17]: multi-objective decomposition with Tchebycheff aggregation, 200 weight vectors, neighborhood size 20, differential evolution operator.

III. Standard MOPSO [20]: baseline MOPSO with adaptive grid archive, without the climate integration, soil model, or adaptive inertia innovations of MOSWARM-CRP.

IV. SPEA2 [33]: strength Pareto Evolutionary Algorithm 2 with fine-grained fitness assignment and density estimation. Population and archive size are 200.

V. Single-Objective GA (SO-GA): a standard GA maximizing only economic return (f1), population size 200, tournament selection, SBX crossover, polynomial mutation.

VI. Rule-Based Rotation (RBR): expert-designed rotation sequences reflecting traditional agronomic best practices for each region. For Mali: millet–millet–fallow–groundnut–sorghum (5-year cycle). For Morocco: soft wheat–fallow–barley–soft wheat (4-year cycle). These represent conventional extension service recommendations.

VII. Random Search (RS): uniform random sampling of feasible rotation plans with constraint repair, serving as a stochastic baseline. 200 × 1000 = 200,000 evaluations (equivalent computational budget).

VIII. All multi-objective methods were configured with equivalent computational budgets: 200 population/swarm size × 1000 iterations = 200,000 fitness evaluations. All methods used the same climate scenario ensemble, crop parameters, and constraint definitions to ensure fair comparison.

Performance was assessed using both standard multi-objective optimization quality indicators and domain-specific agricultural metrics:

Multi-objective quality indicators:

I. Hypervolume (HV): the volume of the objective space dominated by the Pareto front approximation, measured relative to a common reference point. Higher HV indicates better convergence and diversity [60]. The reference point was set to the anti-ideal point across all methods plus 10%.

II. Inverted Generational Distance (IGD): the average minimum distance from each point in a reference Pareto front (constructed as the union of all non-dominated solutions across all methods and runs) to the obtained Pareto front. Lower IGD indicates better convergence and coverage [22].

III. Spread (Δ): a diversity metric measuring the extent and uniformity of the Pareto front distribution [16]. Lower Δ indicates more uniform spread.

Domain-specific agricultural metrics:

I. EAP: mean annual net profit (revenue minus costs) averaged across climate scenarios and rotation years, extracted from the profit-maximizing Pareto front solution.

II. Net Present Value (NPV): discounted cumulative profit over the 10-year planning horizon at a 5% discount rate.

III. Cumulative soil nitrogen depletion (ΔNcum): total nitrogen deficit accumulated over the planning horizon, measured in kg/ha.

IV. Water Use Efficiency (WUE): total crop yield produced per unit of water consumed (kg/m3).

V. Climate Resilience Score (CRS): the worst-case total yield ratio across the K = 10 climate scenarios relative to the historical baseline yield, representing the minimum acceptable performance guarantee.

MOSWARM-CRP parameter settings, determined through preliminary sensitivity experiments on a reduced problem instance, are summarized in Table 4.

HV window for adaptive w τ 20 iterations Penalty coefficient base λ0 10 Independent runs 30

All experiments were conducted on a workstation equipped with an Intel Xeon E5-2680 v4 processor (14 cores, 2.4 GHz) and 64 GB RAM, running Ubuntu 22.04 LTS. The algorithm was implemented in Python 3.10 with NumPy and SciPy for numerical operations. Each independent run (30 per method per case study) was initialized with a distinct random seed.

Table 4 presents the Hypervolume (HV) and IGD indicators across all compared methods for both case studies, reported as mean ± standard deviation over 30 independent runs. Statistically significant differences (Wilcoxon rank-sum test, p < 0.01) between MOSWARM-CRP and each baseline are indicated with superscript markers.

MOSWARM-CRP achieved the highest HV and lowest IGD on both case studies, with statistically significant margins over all baselines. For the Mali case study, MOSWARM-CRP's HV of 0.724 ± 0.018 exceeded NSGA-II (0.651, +11.2%), standard MOPSO (0.667, +8.5%), MOEA/D (0.638, +13.5%), and SPEA2 (0.643, +12.6%). The Morocco case study exhibited a similar pattern, with MOSWARM-CRP achieving HV of 0.698 ± 0.015 compared to 0.623 for NSGA-II (+12.0%). The single-objective GA and RS performed substantially worse on HV and IGD, as expected, since they either collapse or do not explicitly navigate the multi-objective trade-off space. Among the multi-objective baselines, standard MOPSO performed closest to MOSWARM-CRP, suggesting that the PSO search dynamics are inherently well-suited to the rotation planning landscape, and that the domain-specific enhancements in MOSWARM-CRP provide meaningful additional gains.

Table 5 presents the economic performance metrics extracted from the Pareto-optimal solutions. For multi-objective methods, the solution with the highest f1 value on the Pareto front is reported (the "economic-best" knee point). For SO-GA and RBR, the single solution is reported directly.

MOSWARM-CRP delivered the highest Expected Annual Profit (EAP) across both case studies. In Mali, MOSWARM-CRP achieved an EAP of 287.4 × 103 FCFA/ha, exceeding RBR by 23.1% and NSGA-II by 15.1%. In Morocco, MOSWARM-CRP achieved 18.92 × 103 MAD/ha, a gain of 17.9% over RBR and 11.0% over NSGA-II. Notably, MOSWARM-CRP's economic-best solution also outperformed the SO-GA (which maximizes profit exclusively) by 5.9% in Mali and 4.2% in Morocco. This counterintuitive result is attributable to MOSWARM-CRP's soil dynamics model: by maintaining soil fertility through balanced rotations, MOSWARM-CRP sustains higher yields in later planning years, whereas SO-GA solutions tend to deplete soil nitrogen through profit-maximizing but fertility-mining sequences, leading to yield decline over the 10- year horizon.

Table 6 reports environmental performance metrics for the "balanced" Pareto front solution (nearest to the ideal point using normalized Euclidean distance) of each multi-objective method, and the single solution for SO-GA and RBR.

MOSWARM-CRP's balanced Pareto solutions achieved markedly lower cumulative nitrogen depletion and higher WUE than all alternatives. In Mali, ΔNcum was 82.3 kg/ha compared to 142.7 for SO-GA (42.3% reduction), 119.2 for RBR (31.0% reduction), and 97.8 for NSGA-II (15.9% reduction). The pattern in Morocco was analogous, with MOSWARM-CRP reducing ΔNcum by 14.9–39.6% relative to NSGA-II and SO-GA, respectively. These reductions are directly attributable to the embedded soil nutrient dynamics model, which penalizes rotation sequences that create nitrogen deficits and rewards legume inclusion and appropriate fallow scheduling.

WUE was consistently highest for MOSWARM-CRP (Mali: 1.42 kg/m3; Morocco: 2.18 kg/m3), reflecting the algorithm's ability to discover crop sequences and input combinations that extract maximum yield per unit of water consumed. The higher WUE in Morocco reflects the contribution of supplemental irrigation and inherently higher-yielding crop varieties.

Table 4 presents the CRS, defined as the worst-case total yield ratio across the ten climate scenarios relative to the historical baseline. A higher CRS indicates that the rotation plan maintains performance even under the most adverse climate projection.

MOSWARM-CRP achieved the highest CRSs on both case studies. In Mali, MOSWARM-CRP's CRS of 0.812 indicates that even under the worst-case climate scenario (typically the driest SSP5-8.5 projection from CanESM5 or UKESM1-0-LL), the optimized rotation plans maintain 81.2% of historical baseline yield. By contrast, RBRs retain only 60.8% (a 33.6% relative improvement by MOSWARM-CRP), and SO-GA solutions retain only 58.3%. In Morocco, MOSWARM-CRP's CRS of 0.786 compares favorably with 69.5% for NSGA-II (13.1% relative improvement) and 59.2% for RBR (32.8% relative improvement). The explicit minimax resilience objective (f3) and the multi-scenario evaluation framework are the primary drivers of this advantage: by evaluating each candidate rotation across all ten climate futures, MOSWARM-CRP avoids selecting sequences that perform well under average conditions but fail catastrophically under extreme scenarios.

Analysis of the Pareto-optimal rotation plans discovered by MOSWARM-CRP reveals characteristic patterns that differ substantially from traditional practices in both regions.

Mali discovered patterns: The most frequently occurring rotation motifs in the MOSWARM-CRP Pareto archive for the balanced-solution region include: 1) a 4-year cycle of millet → groundnut → sorghum → cowpea, which alternates cereals (Poaceae) with legumes (Fabaceae) in every other year, providing continuous nitrogen fixation inputs and family-based pest/disease break, and 2) a 5-year cycle of millet → groundnut → sorghum → sesame → cowpea, introducing sesame (Pedaliaceae) as a third family for enhanced diversity; and 3) strategic insertion of fallow after cotton in sequences where cotton is included for cash income. These optimized patterns contrast sharply with the traditional millet → millet → fallow cycle common in smallholder practice, which violates same-family succession principles (consecutive millet) and underutilizes the nitrogen fixation potential of legumes. The MOSWARM-CRP patterns achieve higher economic returns through legume diversification (groundnut and cowpea have high market values) while simultaneously building soil nitrogen through fixation, resulting in the synergistic profit–sustainability gains observed in Tables 2 and 3.

Morocco discovered patterns: in the Morocco case study, dominant Pareto-optimal motifs include: 1) durum wheat → chickpea → barley → faba bean (4-year cycle with cereal–legume alternation), 2) soft wheat → faba bean → sunflower → chickpea → barley (5-year cycle with three distinct families), and 3) inclusion of sugar beet every 4–5 years on irrigated parcels for high cash returns, followed by legumes to restore nitrogen. These patterns outperform the traditional wheat → fallow → wheat or wheat → barley → wheat sequences by replacing unproductive fallow years with nitrogen-fixing legumes and by introducing non-cereal crops (sunflower, sugar beet) to diversify both agronomically and economically. The cereals-only traditional rotations, while simple and low-risk, forfeit significant profit from high-value legume and industrial crops and fail to exploit biological nitrogen fixation, resulting in higher fertilizer dependency and greater vulnerability to nitrogen price volatility.

To quantify the individual contributions of MOSWARM-CRP's key innovations, we conducted an ablation study removing each component in isolation. Table 8 reports HV and domain-specific metrics for the Mali case study (results for Morocco follow the same pattern).

The ablation results reveal that each component contributes meaningfully but addresses different performance dimensions. Removing the climate scenario integration causes the largest drop in CRS (from 0.812 to 0.698, a 14.0% reduction), confirming that multi-scenario evaluation is critical for resilience optimization. Without the soil dynamics model, ΔNcum increases by 31.7% (from 82.3 to 108.4 kg/ha), and EAP decreases by 6.6%, demonstrating that accurate soil state tracking is essential for both environmental and long-term economic objectives. Rotations that appear profitable under static soil assumptions lead to yield decline as soil nitrogen is depleted. The adaptive inertia mechanism provides the most modest individual improvement (HV: 0.706 vs. 0.724 for the full model), but contributes to faster convergence and better spread, with particular benefit in the early optimization phases. The combined effect of all three components (full MOSWARM-CRP vs. standard MOPSO) yields a synergistic HV improvement of 8.5%, confirming that the innovations are complementary rather than redundant.

Table 9 reports the mean computational runtime per independent run and the iteration at which 95% of the final HV is achieved (convergence iteration) for each method.

Computational runtime and convergence performance of the compared methods.

MOSWARM-CRP's runtime of approximately 35 minutes is modestly higher than baselines (11.2% more than NSGA-II, 54.9% more than standard MOPSO), primarily due to the multi-scenario fitness evaluation (each particle is evaluated across K = 10 climate scenarios) and the soil dynamics simulation within the fitness function. However, this overhead is modest in absolute terms; 35 minutes is well within the practical tolerance for a strategic planning tool applied once per growing season or planning cycle. Importantly, MOSWARM-CRP converges to 95% of its final HV by iteration 612, substantially faster than NSGA-II (721), MOEA/D (685), and SPEA2 (738), indicating that the adaptive inertia mechanism and PSO social learning dynamics accelerate early convergence despite the higher per-iteration computational cost.