optimLanduse: A package for multiobjective land‐cover composition optimization under uncertainty

How to simultaneously combat biodiversity loss and maintain ecosystem functioning while improving human welfare remains an open question. Optimization approaches have proven helpful in revealing the trade‐offs between multiple functions and goals provided by land‐cover configurations. The R package optimLanduse provides tools for easy and systematic applications of the robust multiobjective land‐cover composition optimization approach of Knoke et al. (2016). The package includes tools to determine the land‐cover composition that best balances the multiple functions a landscape can provide, and tools for understanding and visualizing the reasoning behind these compromises. A tutorial based on a published dataset guides users through the application and highlights possible use‐cases. Illustrating the consequences of alternative ecosystem functions on the theoretically optimal landscape composition provides easily interpretable information for landscape modelling and decision‐making. The package opens the approach of Knoke et al. (2016) to the community of landscape modellers and planners and provides opportunities for straightforward systematic or batch applications.

deciding upon the composition and configuration of land-cover types, given that landscape patterns drive the landscape's ecological value (Duarte et al., 2018) and its ability to satisfy societal demands (Biber et al., 2015;Grass et al., 2020;Kremen & Merenlender, 2018).
Land-cover allocation models have proven helpful in revealing trade-offs between multiple objectives and services provided by different land-cover compositions (e.g. Cavender- Bares et al., 2015).
They allow for the exploration of theoretically optimal land-cover compositions, which satisfy specific sets of requirements on a landscape (Bateman et al., 2013). Optimization procedures, among these particularly Goal Programming and Reference Point approaches, have shown great potential in solving these questions of 'where to put what', (Kaim et al., 2018;Yousefpour et al., 2012), with the landcover shares being the decision variables in the objective function (e.g. Bateman et al., 2013). In landscape planning and environmental modelling, each land-cover type is assigned a set of ecosystem functions, which are predictions that are subject to uncertainties and risks. This is, for example, due to observable natural variability, measuring errors, estimation errors, gaps in knowledge or disturbances (Messerer et al., 2017;Yousefpour et al., 2012).
Incorporating uncertainty into the objective function has a long tradition in landscape modelling. One prominent approach is based on risk reduction through diversification (Markowitz, 1952), which has been applied to landscape portfolios by Flora (1964) and Macmillan (1992), for example. The contributions of land-cover types to ecosystem functions or services are interpreted as asset returns with continuous probability distributions. A risk-reduction effect is achieved if those asset returns are not perfectly correlated.
This has also become an important argument for diversified landscape patterns. A common property of all Markowitz approaches is that they require information about the expectations, uncertainties and also the dependencies between the parameters: that is, any covariances between ecosystem functions across different land-cover types need to be considered. As those parameters are often gathered through Monte Carlo simulations, data requirements in classic Markowitz portfolio optimizations are usually very high, and finding solutions is computationally demanding since only nonlinear solvers are applicable (Knoke et al., 2015).
Recently, robust optimization procedures have become popular for solving portfolio-based land-cover composition problems (Knoke et al., 2015;Yousefpour et al., 2012). Robust approaches do not require any distributional or stochastic information, as the indicator's individual uncertainty spaces are solely defined by their individual best-and worst-case contributions to the land-cover types. The average effect that a land-cover option has on the indicators and one measure of uncertainty are sufficient to construct a range of contribution that a land-cover type may have on an indicator. The optimization program's state space is then constructed by all possible combinations of all best-and worst-case indicator contributions.
Robust solutions guarantee the feasibility of the land allocation over the uncertainty space. The so defined allocation problem is solvable by linear programming, making the solution unique, exact and reproducible. The solver requires no arguments for convergence, and the demands on computational time and resources are very low (Messerer et al., 2017). This formulation also allows for deriving a robust solution that reduces underperformance across a range of multiple indicators. This could be an important attempt towards exploring multifunctional landscapes.
We here present an adaptation of the successfully established and applied robust land-cover optimization approach for multiple objectives by Knoke et al. (2016) for the statistical programming language R (R Core Team, 2022). With this package, we aim to open the method to the community of landscape modellers and planners, and to enhance the possibilities for systematic applications, such as sensitivity analysis, for that approach. We introduce the theoretical optimization framework, sketch the package structure, and explain its functionality using example data from a recently published scientific study.

| ME THODOLOG IC AL BACKG ROUND
The core concept of the robust optimization approach presented here is to find the land-cover composition that minimizes the distance of the worst-performing indicator to its own theoretically best-achievable contribution out of the set of all indicators. We use indicators to describe the relevant functions or services of the landcover types. The best-case contributions and worst-case contributions of an indicator construct individual uncertainty spaces around each indicator and span a broad state space of the optimization program in comparison to approaches without a consideration of uncertainty. The indicators' possible contributions of a land-cover type are bounded by the worst-and best-case outcomes. The robust uncertainty spaces of the indicators are therefore each defined by two parameters. The entire robust uncertainty space of the optimization problem is then framed by all resulting combinations of all indicator contributions. The corners of the state space of the optimization program are determined by all possible combinations of the indicators' worst-and best-case contributions (Knoke et al., 2016).
For each indicator, the number of possible combinations with all indicators' worst-and best-case contributions, also called uncertainty scenarios, is calculated by N U = 2 N L with N L being the number of decision alternatives, that is, land-cover options. For a number of N I indicators, this results in a total number of uncertainty scenarios The solution is the land-cover composition that minimizes the relative distance of the worst-case scenario, for the worst-performing indicator, to its (hypothetical) maximum achievable level. The solution of the program is thus solely determined by the worst-case scenario, that is, the indicator combination that shows under a given land-cover allocation the highest distance to its individually achievable level. Unlike in ordinary goal-programming approaches, the indicators or scenarios with higher performances do not compensate for the worse performing indicators . The program searches the minimum out of the set of maximum distances d iu . d iu are the distances of each indicator to its best-achievable level out of all possible combinations of indicator contributions. Or in other words, the only variable to be minimized is the distance of the worst-performing scenario to its individual maximum achievable level, which we called , thus is not only the objective variable to be minimized, but also determines the maximum threshold of the distances that each indicator in each scenario is allowed to have. This ensures that represents the worstcase scenario. It builds the right-hand side limitations for the constraints of the performances of all indicators , which makes the right-hand side of the optimization program dynamically dependent on the objective.
can be interpreted as the maximum relative distance of the worst-performing scenario to its scenario-individual achievable level out of the set of all scenarios and thus the worst-case scenario of all possible consequences. We thus minimize the worst possible combination of indicators, reflected by the maximum distance. The vector with the scenarios' distances d iu has the length of the number of scenarios N S . It contains the relative differences of each scenario to its theoretically best-achievable level. The distances are calculated as the differences between the maximum achievable level in that scenario, reflected by the most desired indicator contribution for R liu , and the performance of the scenario under a given land-cover composition. Dividing the distances by the range of the lowest and highest uncertainty adjusted indicator min,max iu normalizes the distances.1 As the d iu are calculated by their individual, theoretical best-possible contribution for each indicator i within scenario u as the approach can also be referred to as a Reference Point approach. The R liu are the uncertainty adjusted indicators, which depend on whether the indicator is considered as best case or worst case. As an example, the calculation for indicators where more is better is: with R li being the expectations, and SD li being the uncertainty of the indicators i for the distinct land-cover options l , delivered by the user. f u enables changing uncertainty levels, for example, for sensitivity analysis of the risk attitude of the stakeholders. Figure 3 illustrates the distances d iu for all indicators and uncertainty scenarios (circles) for an explicit example. The maximum of d iu is the relative distance of the worst-performing scenario u of the worst-performing indicator i and thus the worst-case consequence of a land-cover decision under all possible combinations. The difference between the highest and the lowest uncertainty adjusted indicator min,max iu builds the best-achievable level of the scenario iu and, accordingly, can be interpreted as the reference of this scenario. max R liu − R iu , or (1) min .
(2) = max d iu , with i ∈ I, being thesetof indicators and u ∈ U i being thesetof individual uncertainty scenarios for i.
which also means that the min,max iu is calculated with the different uncertainty level f * u . The optimization program can be straightforwardly solved by linear programming. We nevertheless split the solving procedure into a two-stage, nested approach to enable more flexibility for possible further developments. is minimized in the outer optimization (Equation 1). The inner optimization maximizes the distances d iu (Equation 2) regarding the land-cover allocation a l under the known . In its current version, the outer optimization is solved by a self-programmed grid search. The inner optimization is solved by linear programming using the R-package lpSolveAPI (Konis & Schwendinger, 2020).

| OPTIML ANDUSE FUN C TI ON S AND WO RK FLOW
The functions of the package are structured into three successive working steps: initialization, solving and post-processing. The distinct functions network provides higher flexibility and possibly saves calculation time in batch analyses. An initialized object could, for example, be recycled for multiple optimizations. The functions of the third working step contain tools for reappraisal and interpretation of the results (Figure 1).

| Initialization
The initScenario() function integrates the user settings into the data and returns an optimLanduse-object ready for solving. It is capable of processing a long-oriented type of data structure (cf. Gosling et al., 2020;Reith et al., 2020). A detailed example of the data structure required is provided in Chapter 2 in the README of Husmann et al. (2022). The function exampleData() provides an example with the required input structure to be passed as argument Also, a scenarioTable with all scenarios and all calculated interim variables (Equations 3 to 5) is provided.

| Post-processing
The calcPerformance() function calculates and attaches the port-

| E X AMPLE APPLI C ATI ON: AG ROFORE S TRY SYS TEMS
This section illustrates applications of optimLanduse based on the data and results of Gosling et al. (2020). The study took place in a forest frontier region in Eastern Panama and used data from interviews with local farmers. The farmers ranked the performances of different conventional land-cover types and two agroforestry landcover types against various socio-economic and ecological indicators. From the farmer's rankings of each indicator, a mean score and the corresponding standard error of the mean were then calculated for every land-cover type (Table 1).

| Optimizing the land-cover composition with optimLanduse
We applied the optimLanduse functions step-by-step to determine the optimal farm land-cover composition, which minimizes the trade-offs between the different indicators. We imported the example dataset (exampleData()) and initialized an optimLanduse object using initScenario() (Listing 1). This initial object is the basis for the following optimization. We set the uncertainty level (f u , uValue) to 2 and the optimisticRule to 'expectation'. This setting is appropriate to define a risk-averse decision-maker following the philosophy of down-side risk measures; that is, the decision-maker wishes to  . The data in the format required for the optimLanduse package can be accessed via the exampleData( ) function Protecting water a 32 4.0 ± 0.43 4.7 ± 0.41 6.8 ± 0.26 7.6 ± 0.21 7.2 ± 0.44 9.9 ± 0.09 Protecting soil a 32 5.5 ± 0.46 5.0 ± 0.32 6.5 ± 0.40 6.9 ± 0.32 6.6 ± 0.46 9.1 ± 0.38 General preferences a 35 15 ± 3.44 21 ± 3.85 11 ± 3.05 23 ± 3.94 0 ± 0.00 1 ± 0.99 init <-initScenario (coefTable = dat, uValue = 2, optimisticRule = "expectation", fixDistance = NA) Listing 1 Installing the package, loading the example data and defining an initial optimLanduse-object (init). The coefTable passed needs to be strictly formatted in the expected optimLanduse format. More information about this format is given in the help of the initScenario() function and in Chapter 2 of the README .
In the next step, we use the solveScenario() function to solve the initialized optimLanduse object (init) (Listing 2). The resulting land-cover composition is shown in Figure 2. Listing 2 Executing the optimization process and calculating the portfolio's performance.
To analyse how sensitive the optimal land-cover configuration reacts towards increasing risk aversion of the stakeholder, we calculated optimal land-cover portfolios for increasing uValues (Listing 3).
The resulting data frame can be used to compile standard graphics, such as in Figure 4. A more detailed interpretation of the risk attitude is attached in Chapter 4.1 of the README . As a tendency, higher levels of uncertainty lead to more diverse portfolios, since the uncertainty spaces of all indicators increase with increasing uncertainty levels. These broadened individual uncertainty spaces then lead to a broader state space with a higher number of possible candidates for lowest-performing scenarios.

| D ISCUSS I ON AND CON CLUS I ON S
The robust land-cover optimization approach developed by Knoke et al. (2016) is innovative in that it allows for the consideration of uncertainty in ecosystem function provisioning in a multicriteria land-cover allocation model. At the same time, the formulation as a linear programming problem meets the calls for parsimonious models (e.g. Antle et al., 2014). However, so far this approach has been implemented exclusively in Microsoft Excel (e.g. Gosling et al., 2020;Knoke, Kindu, et al., 2020). These implementations have been technically complex, requiring many problem-specific settings and very close-to-data spreadsheet designs, which are time-consuming and prone to errors. As a result, analyses that required model repetitions (e.g. sensitivity analyses) were not feasible. optimLanduse was designed to overcome these limitations, enable fast and systematic parameterizations, and thus make F I G U R E 2 Composition of the optimized farm (based on data shown in Table 1, including all indicators). Each land-cover option is shown in an allocated share (%). optimLanduse provides access to one optimization approach.
Incorporating alternative approaches would, however, certainly be advantageous as it would allow users to choose from a set of approaches or to straightforwardly analyse the differences between these approaches. The package is therefore also meant to serve as a platform for seamless cooperation in the development of landscape modelling approaches. Programmers are invited to copy our code or to fork our project on GitHub as a basis for own landscape-modelling approaches, or to jointly integrate their approaches into optimLanduse.
The robust approach provides one single best land-cover composition under the specific set of preferences out of the entire frontier, using equally weighted decision criteria and uncertainty scenarios. In the current version, all indicators have equal weights. The consideration of differing indicator combinations can still reflect stakeholder preferences. Upcoming versions shall provide the possibility to weigh the indicators so as to reflect the stakeholder's preferences better. These weights could be derived, for example, by an analytical hierarchy process or similar methods as a hybrid multi-criteria approach as, for example, shown by Kaim et al. (2020).
The performances of the indicators give an indication of the indicators to which the result reacts sensitive. See Figure 3 and Chapter 4.2 of the README . The sensitivity towards the risk assumption can be analysed by calculating portfolios with different uValues (f u , Chapter 4.1 of the README).
Necessary parameterizations to calculate the sensitivities have to be conducted by the user. We plan to add a built-in sensitivity report that provides a set of standard sensitivity measures on demand for the indicator expectations, the indicator uncertainties and the uValues. In addition, it may be interesting to analyse how much an indicator can be changed without changing the optimal land-cover composition.
optimLanduse is capable of studying the effects of compositional diversity. With only minor changes, the model could also be used to elaborate the influences of structural diversity on the optimal landcover composition. User-definable restrictions would enable defining groups with different expectations or uncertainty entries of the indicators, which could then be interpreted as spatial regions or any other systematic group.
Our adaptation as an R package has four main advantages: (1) The horizontal and vertical dimensions of the optimization prob-

| 2727
Methods in Ecology and Evoluঞon HUSMANN et al. 1 The package returns a table with all scenarios containing all calculated interim results. See Section 3 for more details.