1. Introduction

Advanced wastewater treatment plants (AWTPs) process wastewater to remove harmful chemicals, nutrients, and microorganisms. By-products of this process are biosolids containing nutrient-rich organic matter, which can then be applied to reuse sites such as local farms, forests, tree farms, and mines for beneficial purposes. The reuse of these bioproducts is carefully regulated by the US Environmental Protection Agency (EPA), to protect humans and their environment from adverse effects (Federal Register, 1993; EPA, 1994a, 1994b).

Recycling of these biosolids products by applying them to the land is just one option. Alternatives include incineration and landfill, but these choices are less desirable from a reuse perspective. Despite some malodorous aspects, land application of biosolids is extensive: 63% in 1998, 66% expected in 2005, and 70% anticipated in 2010 (Oleszkiewicz and Mavinic, 2002), and thus an analysis of how it affects the local population where it is applied is important.

Generally, most of the malodorous nature of biosolids is due to reduced sulphur and nitrogenous compounds (Mosier et al, 1977) such as dimethyl disulfide (DMDS), carbon disulfide (CS₂), and dimethyl sulfide (DMS) (Banwart and Bremner, 1976). These compounds have an odour of rotten cabbage and have low human detection limits (Ruth, 1986). Ammonia (NH₃) and trimethyl amine (TMA) comprise the majority of the nitrogen-based biosolids odours, the former having a medicinal smell and the latter a fishy one (JEA, 1980; Ruth, 1986; Rosenfeld, 1999). Despite the agricultural benefits of properly treated biosolids, there has been some resistance from farm organizations, the food processing industry, as well as local and county groups mostly in agricultural and rural areas. Typical concerns have included the malodorous nature of the biosolids, public health, and the water supply, to name a few areas (De Michele, 2001).

Thus, given the importance and extent of land application of biosolids, proactive planning to mitigate the potential adverse odours is important. Recently, Gabriel et al (2005a, 2005b) have developed several statistical models that predict biosolids odours based on management variables (eg amount of polymer used, number of centrifuges running) and ambient conditions (eg temperature). In Gabriel et al (2005b), the odour scores were measured by field inspectors and provide a surrogate for how the surrounding population would receive these odours. The development of these models was based in part on laboratory experiments that identified key variables relative to biosolids odours (D'Amato and DeHollander, 1999; Rosenfeld, 1999; Sostrand et al, 2000; Kim et al, 2002, 2003; Murthy et al, 2001, 2002; Novak et al, 2002). In addition, in Gabriel et al (2005a), statistical models to predict levels of DMDS were developed based on laboratory experiments.

While statistical models as just described are useful, one needs to also take into account the inherent tradeoffs between resulting biosolids odours and the associated operations and distributional costs. Generally, lower odours will necessitate higher chemical and other costs. Thus, the problem is to produce the least odourous biosolids as inexpensively as possible, while also adhering to certain restrictions (eg do not exceed odour limits at sites). Using the statistical models described in Gabriel et al (2005b), we present a model to jointly minimize odours and costs from the perspective of an AWTP. The resulting multiobjective optimization model can then be used at AWTPs to balance these competing objectives and provide quantitative measures of the tradeoffs. Such a model can estimate, for example, if a reduction of biosolids odours by 10% is desired, how much the optimal costs would have to increase to make this happen. This approach takes into account stakeholders in the decision-making process as advocated in (De Michele, 2001) in that the odour effects on the surrounding population as well as the plant and distribution costs are considered. It necessarily goes beyond work that just considers one objective such as minimizing the total costs for land application of biosolids (Crohn and Thomas, 1998).

Many studies and published papers have dealt with environmental and/or waste treatment issues. For relevance to the current paper, we focus only on a selected subset of studies that have used either mixed integer or multiobjective programming. An early work by Escudero (1978) provides a comparative analysis of methods for representing nonlinear functions of a single argument in mathematical programming for pollutant emission control policies and standards. Ciric and Huchette (1993) applied multiobjective programming techniques to analyse the sensitivity of maximum net profits of a chemical process to changes in the waste treatment cost when discrete terms need to be optimized. Later, Chang et al (1997) investigated optimal wastewater treatment strategies for water pollution control in a river basin using interactive, fuzzy, interval, multiobjective, mixed integer programming. Another related study by Chang and Wang (1997) considered a fuzzy goal programming approach for optimal planning of solid waste management systems in a metropolitan region. Papalexandri and Dimkou (1998) proposed a decomposition algorithm to identify the noninferior solution set of a multiobjective, mixed integer, nonlinear program through a parametric optimization problem. The proposed algorithm was applied to the process synthesis and simultaneous product/process design problems. Later, Shih and Cheng (2001) studied hauling costs and environmental impacts when considering two competing objectives in transportation planning for industrial wastes. Nema and Gupta (2003) presented a multiobjective integer goal programming model for analysing problems involving planning and design of regional hazardous waste management systems. Lastly, Marseguerra et al (2004) used a multiobjective genetic algorithm approach to model the optimization of technical specifications of a nuclear safety system.

The current work distinguishes itself from the previous studies in that a detailed model of both the wastewater treatment processes as well as distributional aspects are considered under the framework of a multiobjective optimization. In this optimization, the objectives of biosolids odour levels and total wastewater treatment and biosolids distribution costs are balanced to arrive at a set of Pareto-optimal decisions. As far as we know, this simultaneous treatment of multiobjective optimization of wastewater and biosolids management is novel. Thus, it represents an extension to the environmental modelling efforts described above.

In addition, the current paper distinguishes itself in the way that it handles certain bilinear, that is, nonconvex, functions relating to the product of chemical levels and processing variables. Based on Schur's decomposition (Horn and Johnson, 1985) and the use of special ordered set (SOS2) variables (Beale and Tomlin, 1970; Beale and Forrest, 1976; Escudero, 1978), the problem is closely approximated by a mixed integer linear program. See Gabriel et al (2006) for a discussion of the practical advantages of this approach as compared to other methods.

The rest of this paper is organized as follows. In Section 2 we describe the wastewater treatment processes at the Blue Plains AWTP used as a case study to validate and illustrate the proposed model. In Section 3 we describe a multiobjective optimization model to jointly minimize biosolids odours and costs and in Section 4 present some representative results using the Blue Plains AWTP database. In Section 5, we summarize the contributions of this paper and describe possible future work. Lastly, the appendices contain relevant mathematical details.

2. Biosolids processing and statistical analysis

The Blue Plains AWTP in Washington, DC, operated by the DC Water and Sewer Authority (DCWASA) serves Maryland, Virginia, and the District of Columbia with approximately 40, 10, and 50% of the processing, respectively in these areas. DCWASA receives, on the average, 370 MGD (1402.3 MLD) of wastewater/storm flow through a series of preliminary, primary, and secondary processes with approximately 1200 wet tons of biosolids produced daily. These biosolids are then distributed to some 3000 fields (eg farms) in Maryland and Virginia.

This flow initially passes through bar screens for trash removal and is then treated with iron salts (FeCl₃) for phosphorus removal. The wastewater continues to flow to the primary sedimentation tanks.

Primary process: Organic suspended solids and phosphorus are removed from the wastewater flow by slowing the flow and allowing gravity settling in primary sedimentation tanks. Gravity thickened solids are pumped to the blend tank. The primary wastewater effluent flows to secondary reactors.

Secondary process: Primary effluent flows to the secondary aeration reactors and is mixed with waste pickle liquor FeCl₃ for additional phosphorus removal. The amount of pickle liquor added is one of the independent variables that could affect downstream biosolids odours and was considered in the statistical modelling. This mixture (mixed liquor) is aerated in the secondary reactors and aerobic microorganisms are grown at a high rate to remove suspended and colloidal carbon, and phosphorus. The mixed liquor from the aeration reactors flows to the secondary sedimentation tanks and the solids are gravity settled. A portion of these settled solids (RAS) are returned to the aeration reactor to maintain microbial proportions. The rest of the settled solids are pumped to dissolved air floatation thickening tanks. The secondary effluent flows to the nitrification/denitrification process.

Nitrification/denitrification: The nitrification/denitrification process removes ammonia from the secondary effluent using methanol as a food source for aerobic and anoxic microbe to remove N₂ gas. The settled solids are pumped to either the reactors or to dissolved air floatation thickening tanks.

Dissolved air floatation (DAF): Settled solids from the secondary and nitrification/denitrification sedimentation tanks are pumped to air floatation tanks and mixed with compressed air and polymer to coagulate and thicken the solids. After thickening the solids are pumped to the blend tank.

Blend tank: Gravity thickened primary solids and DAF thickened secondary are mixed together in the blend tank.

Dewatering and lime stabilization: The solids from the blend tank are mixed with polymer and then dewatered by high solid centrifuges and belt presses. The dewatered cake is then mixed with lime (CaO) for pathogen reduction to EPA Class B level.

A more detailed description of these processes and potential independent variables used to predict odours via statistical models is given in (Gabriel et al, 2005b); an example of one of these statistical models is given below.

where Y_i is the Inspector C's biosolids odour score on day i in the winter period, and X_j is the jth independent variable with the independent variables in order as:

• The minimum temperature on day d-1 (F)
• Blanket depth at the secondary east tank on day d-1 (feet)
• The number of contractor belt filter presses in service on day d-1
• The number of contractor centrifuges in service on day d-1
• The amount of lime used on day d-1 (lbs/dry ton)
• A dummy variable for when the polymer used on day d-1 was greater than 200.05 lbs/dry ton
• A dummy variable for when the lime used on day d-1 was less than 308 lbs/dry ton

The notation 'd-1' means 1 day before the day the odour score was recorded. Equation (1) contains important variables that can be used by AWTP management. For instance, when the number of centrifuges operated by the contractor (as opposed to DCWASA) is increased by one, all things being equal, the odour score should drop by about 1.91 points (presumably due to more efficient processing). Thus, management should ensure that as many of these centrifuges are running as possible. Clearly, there is a cost associated with such a management decision. The multiobjective optimization model described below allows one to make a quantifiable and comprehensive tradeoff between odour levels and costs by varying such things as the number of DCWASA centrifuges. Of course, other factors such as the ambient temperature are not under control of the AWTP but are important for calibrating the model and developing the tradeoffs between biosolids odours and costs.

3. Multiobjective optimization model: case study for Blue Plains AWTP

The multiobjective optimization model to minimize biosolids odours as well as costs operates for a given time horizon in which each day is denoted by dD. For example, if one week is chosen, D={1, 2, 3, 4, 5, 6, 7}. During this timeframe, the biosolids will be distributed to a set of reuse fields F with each field denoted as f. All of the 3000 or so fields under DCWASA's jurisdiction could be used but typically only a subset of these will be included for computational as well as management reasons. The set of distribution contractors hauling the biosolids to these fields is given by C with each contractor denoted as cC. Conversely, there is also one on-site processing contractor responsible for a portion of the dewatering process at the Blue Plains AWTP. Variables and data indexed by k refer to this on-site contractor. This contractor uses both belt filter presses and centrifuges to dewater the sludge, whereas the remaining amount of material is dewatered by DC WASA using just centrifuges.

The main variables in the model are

F_d^k

Percent flow from blend tank to the on-site contractor on day d,

F_d^dc

Percent flow from blend tank to DCWASA on day d,

F_dcf

Amount of biosolids applied to field f by contractor c on day d (tons),

L_d^dc

The amount of DCWASA lime used on day d (lbs/dry ton of biosolids),

LD_d

A dummy variable for when the lime amount on day d is less than 308 lbs/dry ton biosolids (0 or 1); note that the 308 lbs/dry ton represents a 20% fractile value,

C_d^k

The number of contractor centrifuges in service on day d,

B_d^e

Blanket depth for secondary east tank on day d (feet); note that the blanket depth is the level of solids that build up in the sedimentation tanks in the secondary process,

W_d^dc

Wet tons of biosolids produced by DCWASA on day d,

W_d^k

Wet tons of biosolids produced by on-site contractor on day d,

O_d

Odour of biosolids on day d.

The last four variables are intermediate variables computed from the others in the following way.

Note that the relationship between blanket depth at the secondary east tank at Blue Plains and the number of contractor centrifuges in service in (2) was determined from year 2002 data. This relationship indicates the greater the number of centrifuges running the lower the corresponding odour. Also, the last day of time horizon in question (|D|) was not included here because biosolids' odour on each day was determined from the blanket depth the day before. In other words, the blanket depth on the last day would not affect the optimal solutions in the model.

Wet tons of biosolids W_d^dc and W_d^k in Equations (3) and (4) were determined from the percents of solids in dry-tons biosolids S^dc and S^k, respectively.

Equation (5) corresponds to (1) but we changed the variable notations to correspond to those used in the optimization model.

The main data used in the model were

B_d

Dry tons of biosolids on day d (tons),

B_d^daf

Dry tons of biosolids in the DAF tanks (tons),

G_d^k

The number of contractor belt filter presses in service on day d,

S^kf

Percent solids in biosolids produced by the on-site contractor,

S^dc

Percent solids in biosolids produced by DCWASA,

H_cf

Hauling cost to field f by contractor c ($/ton),

F_f^up

Field capacity (tons),

P_d^daf

The amount of polymer DCWASA adds in the DAF tanks on day d (lbs/dry ton biosolids),

P_d^de

The amount of polymer DCWASA adds in the dewatering process on day d (lbs/dry ton biosolids),

PD_d

A dummy variable for when the sum of DAF polymer and DCWASA dewater polymer used on day d was greater than 200.05 lbs/dry ton biosolids (0 or 1); note that the 200.05 lbs/dry ton represents the 80% fractile value,

T_d^min

Minimum temperature on day d (degrees Fahrenheit).

The two objectives to be minimized simultaneously are the total odour

and the total processing and distribution costs

where PC_d^dc is the polymer cost on day d ($/day), with

PC=Polymer cost ($/lb)

LC=Lime cost ($/lb)

CC_d^k=Contractor centrifuge's maintenance and operating cost on day d ($/day), with

CC=Cost (maintenance cost and operating cost) of keeping one centrifuge running ($/centrifuge); note that this was a cost that was varied in the sensitivity analysis

K_d^O=on-site contractor operating cost on day d with

where F_d^kB_d represents dry tons of biosolids processed by the on-site contractor on day d. The on-site contractor operating cost was $71.75 per dry ton of biosolids for the first 150 dry tons and then $50.25 per dry ton of biosolids from 150 to 250 tons as described by DCWASA management. This was the rationale for the step function (11).

The on-site contractor lime cost was calculated from pre-lime and post-lime, where pre-lime was $ 10.1 per dry ton biosolids and post-lime cost was $7.4 per dry ton biosolids. Note, however, that only 67% of biosolids processed by onsite-contractor was pre-limed. Adding these two lime costs, we obtained (12).

The decision variables are subject to certain constraints. First, there is conservation of flow (%) between DCWASA and the on-site contractor as well as nonnegativity of the flow percentages:

Since no storage of biosolids at the plant is modelled, conservation of biosolids means that

This equation can be explained as follows. First, L_d-1^dcB_d-1/2000 represents the additional wet tons of biosolids due to lime addition; B_d-1^dafP_d-1^daf/2000 and B_d-1P_d-1^deF_d-1^dc/2000 represent the additional wet tons of biosolids due to polymer addition at the DAF and dewatering processes, respectively. Therefore, the sum of the first five terms in (15) represents total hauling weight.

The amount of biosolids distributed over the time horizon to each field cannot exceed the given capacity (see Table 1). These restrictions as well as the nonnegativity of F_dcf are given as follows:

Table 1 - Sample of reuse field data.

Full table

The daily amount of biosolids hauled to reuse fields by each contractor c is limited by CL_c^up with the constraint given as follows:

The next set of constraints enforces additional logic to show that for each reuse field fF, if any distribution is made, that is, F_dcf>0, then the odour limit for that field should not be exceeded.

where O_f^up is the upper bound on odour level for each f (odour level).

To see the logic of these three constraints, note that when F_dcf>0, y_df must equal one in (19) implying that O_dO_f^up in (20). Also, when y_df=1, the bound of 1200 relates to the maximum daily amount of biosolids that can be distributed by a contractor to a field. When F_dcf=0, the binary variable y_df can equal zero or one. When y_df equals zero, no shipment is made and the odour score is not relevant. Moreover, no realistic additional restrictions are imposed on the odour variable O_d in (20) in this case as well since the value of 1200 is way beyond any valid odour score.

Additionally, if F_dcf>0, it should be in the range [30, 1200] to reflect actual distribution values. The following two constraints are meant to increase the effectiveness of transportation of biosolids from DCWASA based on the capacity for each truck being approximately 30 tons. As stated above, 1200 tons represent the maximum amount of each shipment.

In addition, there are minimum and maximum values for the amount of lime used and the number of centrifuges that can be running, given respectively as

where L^low is the lower bound on the amount of lime DCWASA added (lbs/dry ton biosolids), L^up is the upper bound on the amount of lime DCWASA added (lbs/dry ton biosolids), C^kup is the upper bound on the number of contractor centrifuges running, C^klow is the lower bound on the number of contractor centrifuges running.

The integrality constraints for the number of contractor centrifuges running is also imposed.

Next, the percentage flow from the blend tank to the on-site contractor on day d should correspond to the number of contractor's centrifuges and the number of contractor's belt filter presses in service. The next constraint makes this correspondence between the sum of the total number of DCWASA's centrifuges, and the total number of contractor's centrifuges and belt filter presses.

Thus, the full multiobjective optimization model is as follows:

s.t. (2)–(27).

Problem (28) as stated has a linear objective function (f₁) as well as a bilinear one (f₂). The bilinearity for the second objective is a result of the product of the two decision variables L_d^dc and F_d^dc in the term LC_d^dc=L_d^dcF_d^dcB_dLC,dD{|D|}. Bilinear functions are not convex and make finding a solution challenging due to this lack of convexity. These bilinear terms in the second objective function are approximated via Schur's decomposition (Horn and Johnson, 1985) and the use of SOS variables of type 2 (SOS2) (Beale and Tomlin, 1970; Beale and Forrest, 1976) using a technique described in (Gabriel et al, 2006) resulting in linear and binary restrictions. Appendix A shows mathematical details on Schur's decomposition and the SOS type 2 variables used. Also, Appendix B demonstrates how we decomposed the bilinear term in our problem explicitly. All the other constraints described above are linear with the exception of binary restrictions on some of the variables. Thus, the resulting problem is still nonconvex but is a mixed integer linear program for which the constraint method can be applied to generate Pareto-optimal points.

4. Case study for the Blue Plains Advanced Wastewater Treatment Plant

A case study was performed to analyse the operations at DCWASA the world's largest advanced wastewater treatment plant serving the Washington, DC metro area. For the case study, the time horizon was 3 days using 2 January 2002–4 January 2002 as a typical time period. Additionally, the model considered three hauling contractors and 40 reuse fields. The data from these fields were provided by Maryland Environmental Services (MES). Selected data for eight of the fields including their tonnage biosolids capacity at the end of August 2004 are displayed in Table 1.

4.1. Population and school data

To aid in the decision of where to send the biosolids, data on the local population and number of adjacent schools were collected. All things being equal, the odour thresholds for reuse sites should be lower for more densely populated areas and/or locations with close proximity to schools. Population densities (people/square mile) were obtained from Environmental Systems Research Institute, Inc. (ESRI)'s Data & Maps CD-ROM. ESRI's Data & Maps CD-ROM provides population densities at county, track, and block levels. Tracks are small, relatively permanent statistical subdivisions of a county or a statistically equivalent entity in accordance with US Census Bureau guidelines. Blocks are bounded on all sides by visible features such as streets, streams, and railroad tracks, and by invisible boundaries such as city, town, and county limits.

In order to assign a population density to each reuse field, each field's location was plotted on the Maryland and Virginia block population density map layer based on its latitude and longitude (see for example Figure 1 for Louisa County, Virginia). Note that each block may contain multiple reuse fields. Then the population density of each block was assigned to each reuse field contained in that block.

Figure 1.

Reuse fields, schools and block population density, Louisa County, Virginia.

Full figure and legend (217K)

Locations of all schools in Maryland and Virginia were also obtained from ESRI's data and map CD-ROM. Table 2 shows the block population density of the eight representative fields and the number of schools in a radius of 3 miles, where OID number refers to the same OID as in Table 1.

Table 2 - Reuse fields with block population density and number of schools within 3 miles.

Full table

4.2. Temperature data

Temperature data were collected from the National Climate Data Center (NCDC) website (http://www.ncdc.noaa.gov/oa/ncdc.html) consistent with Gabriel et al (2005b). To represent the temperature at the Blue Plains AWTP, data on the temperature at the National Airport Station, data for the closest station of NCDC were collected. Three types of temperature data were used: minimum, maximum and average values. Temperature data from 2 January 2002 to 4 January 2002 shown in Table 3 were used as input data for the optimization model.

Table 3 - Temperature on selected days (degrees F).

Full table

4.3. Blanket depth data

The depth in feet, of the suspended solids in the secondary settling tank, called the blanket depth was identified in Gabriel et al (2005b) as an influential variable relative to biosolids odours. This variable roughly determines the capacity of waste activated sludge in the tanks. At the Blue Plains AWTP, there are three types of blanket depth: blanket depth secondary west even tanks, blanket depth secondary west odd tanks, and blanket depth secondary east tanks. The blanket depth data from 1 January 2002 to 8 January 2002, that were used as representative input to the optimization model are shown in Table 4.

Table 4 - Blanket depth in feet on selected days.

Full table

4.4. Operating and maintenance cost data

Data related to the operations and maintenance of the centrifuges were very difficult to obtain. Since the number of operating centrifuges can be an important factor in determining the biosolids odour levels, these variables were crucial to the optimization model. Consequently, these costs were treated as a parameter to be varied to see their effect on the overall results. A cost of $196 per day to operate and keep one centrifuge running was assumed as the nominal value based on conversations with DCWASA management.

4.5. Odour thresholds for each reuse field

The biosolids odour scores were determined subjectively by field inspectors with the range of scores from 0 (least malodorous) to 9 (most malodorous). The odour thresholds for each reuse field f were not available from the inspectors and were thus computed from scratch. In particular, many choices for these values are possible. However, since the malodorous aspects generally affect people, factors that involved the local population should be used. Also, all things being equal, the more people potentially affected by the odours (within a given radius), the lower the threshold should be.

Using these guiding principles, each odour threshold value O_f^up was computed as follows. The inputs for this odour threshold index calculation were the number of schools in a radius of 3 miles and the population density for the block where the field was located as shown in Table 2. The steps to calculate the odour threshold index were as follows:

(1) For a field that had no schools in a radius of 3 miles, assign a school index of 9 to that field (relating to the maximum possible odour score given by the inspectors).

(2) For a field that had a nonzero number of schools in a radius of 3 miles:

(a) Calculate School_Score1 as

with the coefficients 2.5, 1.5, and 0.5 representing the relative importance of the associated number of schools in each group. Larger values for School_Score1 indicate a higher value of schools in proximity to the reuse field that might be affected by the odours. Taking the reciprocal of this score would relate to the threshold value as shown next.

(b) Calculate School_Score2 as

The values for School_Score2 needed to be normalized so that they went from 0 (lowest threshold) to 8 (highest threshold apart from Step 1); this is done in the next step.

(c) Calculate School_Index as

A similar line of reasoning was applied to the population density values for the area near the reuse field. In particular, the following computational steps were used:

(1) Calculate the population index as

(2) Normalize this value so that the value is between 0 and 9

The last step was to combine the school and population indices by averaging them.

4.6. Computational results

The optimization model was solved using the XPRESS-MP optimization solver (version 15.25.02) with the modelling interface XPRESS-IVE (version 1.15.04) and the MOSEL language (version 1.4.1) (www.dashoptimization.com). The computer used was a Dell OPTIPLEX GX270, Pentium 4 with a CPU of 3.00 GHz and 2 GB of RAM memory. Three cases were examined: 'Base', '$296', and 'KCost1'. Solving problem (28) constituted the Base case. The other two cases were distinguished from the Base Case by varying the contractor centrifuge cost and the contractor piecewise linear operation costs. In particular, the '$296' case used $296 instead of an cost of $196 in the Base case for the operations and maintenance of the contractor centrifuge. The 'KCost1' case differed from the Base case in that the contractor's piecewise linear operations cost function was changed to

instead of function (11) resulting in a different and higher cost structure. The motivation for using these three cases was to see how the tradeoffs between odour reduction and costs changed by varying the parameters specified above.

In generating the Pareto-optimal set, we tried the weighting and constraint methods (Cohon, 1978). Although the weighting method is simple to implement, for this case study, only one Pareto-optimal point was obtained out of 38 subproblems. In fact, nine subproblems were solved to optimality but all of them yielded the same Pareto-optimal point. Another interesting fact was that all nine Pareto points were found within some milliseconds. Moreover, the solver could not finish the branch and bound search within a preset maximum running time for the other 29 subproblems (the preset maximum running times were arbitrarily chosen and were 2 h for one subproblem and 1 h and half an hour for some subproblems). This was undoubtedly due to the nonconvex nature of the models and, hence the existence of 'duality gap' points. The weights w₁, optimization status, running time, total cost, and total odour for all 38 subproblems are shown in Table 5.

Table 5 - Optimal solutions for the Base Case using the weighting method.

Full table

The constraint method faired better in finding Pareto-optimal solutions presumably due to its superior ability to find duality gap points, however, with a large computational burden (ReVelle and McGarity, 1997). Using this method, several subproblems were solved varying the maximum value . First, the total cost was minimized without constraining the total odour and the resulting total odour was recorded. Using this value as the maximum value , we next solved the first subproblem by minimizing the total cost. Then, the next subproblems were solved by decreasing the total odour constraint right-hand side by 0.1 until the problem was infeasible. One of the criteria to guarantee that the optimal solution obtained was Pareto optimal was that the constraint on the total odour should be binding (Cohon, 1978). In addition, it turned out that for some subproblems, the solver could not find the solutions within the preset maximum times. Again, this was due to the nonconvex nature of the models as we also encountered in the weighting method. For these subproblems, we recorded the total costs yielded when the search stopped and used each recorded total cost as the maximum right-hand side value to generate each additional subproblem with minimizing total odour as the objective. Using this approach, we were able to find several more Pareto points.

The solution times for finding the Pareto optimal points varied. Sometimes solutions were obtained in milliseconds while significantly longer times were needed in other cases. Examples of specific solution times for selected runs can be found in Tables 6, 7 and 8, for, the Base, $296, and KCost1 cases, respectively. In addition, note that the long run time of 1767 s in Table 6 is due to the nonconvex nature of the problem.

Table 6 - Running time for generating Pareto-optimal points for the Base case using the constraint method.

Full table

Table 7 - Running time for generating Pareto-optimal points for the $296 case using the constraint method.

Full table

Table 8 - Running time for generating Pareto-optimal points for the KCost1 case using the constraint method.

Full table

4.7. Analysis of the three cases

Figure 2 shows three sets of Pareto-optimal points generated from the Base, '$296', and 'KCost1' cases. First, note that from this figure, certain odour scores were negative. This was due to the fact that the data in the time horizon selected went beyond the range of data used to generate the statistical model for odour prediction. However, the odour scores were only relative values and therefore as stated, are useful for management decisions.

Figure 2.

Approximation of the Pareto-optimal set for the three cases.

Full figure and legend (52K)

Increasing the centrifuge costs by $100 had relatively little effect on changing the approximation to the Pareto-optimal set of decisions as compared to the Base case. In comparing Pareto-optimal points common to both the Base and '$296' cases that had the same odour level, apparently the values of the optimal decision variables were the same. Consequently, the differences in cost were attributable to just the centrifuges. The extra $100 in the '$296' case increased the total costs by that amount times the number of centrifuges chosen to be operated. As can be seen in the upper left portion of Figure 2, when there were high odour levels, there were no centrifuges selected, hence the Pareto-optimal points were the same for these two cases.

The effect of changing the contractor operations cost in the 'KCost1' case was similar to what was described above. In particular, the costs of polymer, DCWASA lime, running the on-site contractor centrifuges, on-site contractor lime, and distribution, were apparently the same as in the Base case, given that the Pareto points with the common odour were compared. The difference in costs is thus attributable to the increased contractor operations costs shifting the Pareto curve up. Thus, for common odour values, it can be seen that changing either the contractor centrifuge or contractor operations costs did not affect the optimal decision variable values.

Another interesting way to analyse these three cases is the estimate of the number of dollars needed to reduce the odour by one index point. To obtain these numbers, statistical regressions on the Pareto points of these three cases were done. More specifically, these regression results were as follows:

Each of these regressions had R² and adjusted R² values >0.99 with statistically significant coefficients (ie, t-statistics greater than 85 in absolute value). From these equations, to reduce the odour by one index point, one needed to pay on average 1/0.0011=$909, 1/0.0011=$909, and 1/0.0009=$1,111, respectively, for the Base, '$296' and 'KCost1' cases. Figures 3, 4 and 5 show three sets of Pareto-optimal points and their corresponding regressions generated from the Base, '$296', and 'KCost1' cases, respectively.

Figure 3.

Regression for the approximation of the Pareto-optimal set, Base case.

Full figure and legend (50K)

Figure 4.

Regression for the approximation of the Pareto-optimal set, '$296' case.

Full figure and legend (46K)

Figure 5.

Regression for the approximation of the Pareto-optimal set, 'KCost1' case.

Full figure and legend (47K)

In trying to understand the cost-odour relationship better, it is important to delve a bit more deeply into the analysis. To this end, we consider just the Base and '$296' case in the next analysis. What is important are the ratios of the decrease in odour due to spending one dollar on either lime or centrifuges. The ratio for the decrease in odour due to dollars spent on lime was 0.01/0.06=0.17 for both the Base and '$296' cases, obtained by taking the coefficient for lime in the statistical model (0.01) and dividing it by the dollars/pound of lime (0.06). Similar ratios for centrifuge dollars were 1.91/196=9.75E^-3 and 1.91/296=6.45E^-3 for the Base case and '$296' cases, respectively.

In considering just these two cases, there were three Pareto points with the same odour and cost values: (7.1824, $222241.5935), (7.0825, $222412.1594), and (6.8825, $222848.9660). In spite of the differences in centrifuge costs between these two cases, the reason these three points were common to both can be explained as follows. First, when we minimize the cost and constrain the total odour objective, more odour would be decreased per dollar spent on lime than centrifuges, so it is more efficient to use as much lime as possible. In fact, the model met the odour objective constraint while using less than the maximum allowable amount of lime and no contractor centrifuges were needed. Next, when we minimize odour while constraining cost it is more efficient to spend the limited budget on a resource that could reduce odour more, that is, lime. Again, no contractor centrifuges were needed as reported by the model. Thus, these points were common to both cases since no centrifuges were selected and these cases only differed in centrifuge costs. However, these conclusions change for points with lower total odour for which centrifuges were needed.

The last item to be analysed was the biosolids distribution pattern. As reported by the model, on any given day, biosolids were assigned to contractor 3 until its daily capacity limit of 600 tons per day was reached. Then the rest of the biosolids were assigned to contractor 1. This distribution pattern resulted since the hauling cost for contractor 3 was cheapest among all three contractors. Also, contractor 2 was not used at all because it had the most expensive hauling costs.

As for the distribution to reuse fields, biosolids were sent to any field where its odour threshold limit was higher than the computed odour score for that day and for which the field capacity was higher than 30 tons. These 30 tons were affected by the distribution constraint that any F_dcf must be at least this amount. In addition, multiple distribution solutions were encountered by first solving the model, then constraining a non-zero F_dcf in the solution to zero and then rerunning the model. Multiple solutions can be explained as follows. First, the distribution cost of the biosolids was independent of where they were sent (this is true in practice). More specifically, the distribution cost was charged on a per tonnage basis. Thus, there were multiple choices of where to send the biosolids for the same cost.

5. Summary and future work

In this paper, we have presented a novel multiobjective optimization formulation for managing both biosolids odours and costs at advanced wastewater treatment plants (AWTPs). It is anticipated that AWTP management could make use of this model to proactively handle odour complaints at reuse sites close to populations while at the same time being cost efficient about their operations. We have validated the model and run it on a set of cases to demonstrate the tradeoffs between odour levels and the resulting costs. It is fair to say that such an analysis is rarely done, if at all, at AWTPs given the complexity of the task. However, the results of this sort of analysis are quite useful for management decisions.

We anticipate that future work will concentrate in two directions. First, algorithms to more efficiently solve the subproblems that arise from using the constraint method (eg Dantzig–Wolfe decomposition (Dantzig, 1963)). This will allow for generating more Pareto-optimal points as well as for solving larger problems (eg longer planning horizons). Second, we anticipate that additional realistic features will be added (eg taking into account forecasted wind speed at reuse sites) and that a geographic information systems interface will also be developed.

References

1. Banwart WL and Bremner JM (1976). Evolution of volatile sulfur compounds from soils treated with sulfur-containing organic materials. Soil Biol Biochem 8: 439–443. | Article | ChemPort |
2. Beale EML and Forrest JJH (1976). Global optimization using special ordered sets. Math Prog 10: 52–69. | Article |
3. Beale EML and Tomlin JA (1970). Special facilities in a general mathematical programming system for non-convex problems using ordered sets of variables. In: Lawrence J (ed). Proceedings of the Fifth IFORS Conference. Tavistock, London, pp 447–454.
4. Chang NB, Chen HW, Shaw DG and Yang CH (1997). Water pollution control in river basin by interactive fuzzy interval multiobjective programming. J Envir Eng 123: 1208–1216. | Article | ChemPort |
5. Chang NB and Wang SF (1997). A fuzzy goal programming approach for the optimal planning of metropolitan solid waste management systems. Eur J Opl Res 99: 303–321. | Article |
6. Ciric AR and Huchette SG (1993). Multiobjective optimization approach to sensitivity analysis: waste treatment costs in discrete process synthesis and optimization problems. Ind Eng Chem Res 32: 2636–2646. | Article | ChemPort |
7. Cohon JL (1978). Multiobjective Programming and Planning. Academic Press: New York.
8. Crohn DM and Thomas AC (1998). Mixed-integer programming approach for designing land application systems at a regional scale. J Environ Eng ASCE 124: 170–177. | Article | ChemPort |
9. D'Amato II RM and DeHollander GR (1999). Gaseous emissions from wastewater facilities. Water Environ Research 71: 715–720. | ChemPort |
10. Dantzig GB (1963). Linear Programming and Extensions. Princeton University Press: Princeton.
11. De Michele E (2001). The national biosolids partnership: background, progress, and the future. In: Proceedings of Biosolids Management in the 21st Century. Water Environmental Federation, Alexandria, VA (CD-ROM).
12. Environmental Protection Agency (EPA) (1994a). Biosolids recycling: beneficial technology for a better environment. Report no. EPA 832/R-94-009, U.S. Environmental Protection Agency, Office of Water, Washington, DC.
13. Environmental Protection Agency (EPA) (1994b). Plain English guide to the EPA part 503 biosolids rule. Report no. EPA/832/R-93/003, U.S. Environmental Protection Agency, Office of Wastewater Management, Washington, DC.
14. Escudero L (1978). A comparative analysis of linear fitting for non-linear functions on optimization, a case study: air pollution problems. Eur J Opl Res 2: 398–408. | Article |
15. Federal Register (1993). Standards for the use or disposal of sewage sludge: final rules. Federal Register 58: 9248–9415.
16. Gabriel SA et al (2005a). Prediction of dimethyl disulfide levels from biosolids using statistical modeling. J Environ Sci Health, Part A 40: 2009–2025. | ChemPort |
17. Gabriel SA, Vilalai S, Peot C and Ramirez M (2005b). Statistical modeling to forecast odor levels of biosolids applied to reuse sites. ASCE J Environ Eng (in press).
18. Gabriel SA, García-Bertrand R, Sahakij P and Conejo A (2006). A practical approach in approximating bilinear functions in mathematical programming problems by using Schur's decomposition and SOS type 2 variables. JORS (doi:10.1057/palgrave.jors.2062052). | Article |
19. Horn RA and Johnson CR (1985). Matrix Analysis. Cambridge University Press: Cambridge.
20. Japanese Environmental Agency (JEA) (1980). Reports of studies on the measurement of offensive odors (from 1972–1980) (in Japanese), Code D3, Tokyo.
21. Kim H et al (2002). Characterization of wastewater and solids odors using solid phase microextraction at a large wastewater treatment plant. Water Sci and Tech 46(10): 9–16. | ChemPort |
22. Kim H et al (2003). Examination of mechanisms for odor compound generation during lime stabilization. Water Environment Research 75: 121–125. | Article | PubMed | ChemPort |
23. Marseguerra M, Zio E and Podofillini L (2004). A multiobjective genetic algorithm approach to the optimization of the technical specifications of a nuclear safety system. Rel Eng and Syst Safety 84: 87–99. | Article |
24. Mosier AR, Morrison SM and Elmond GK (1977). Odors and emissions from organic wastes. In: Elliott LF and Stevenson FJ (eds). Soil for Management of Organic Waste and Waste Waters. Soil Science Society of America, Madison, Wisconsin, pp 532–569.
25. Murthy SN et al (2001). Mitigation of odors from lime stabilized biosolids. WEF Residuals and Biosolids Management Conference, San Diego, CA (CD-ROM).
26. Murthy S et al (2002). Mechanisms for odour generation during lime stabilization. IWA Biennial Conference, Melbourne, Australia (CD-ROM).
27. Nema AK and Gupta SK (2003). Multiobjective risk analysis and optimization of regional hazardous waste management system. Pract Periodical of Hazard, Toxic, and Radioactive Waste Mgmt 7: 69–77. | Article | ChemPort |
28. Novak J, Glindemann D, Murthy SN, Gerwin SC and Peot C (2002). Mechanisms for generation and control of trimethyl amine and dimethyl disulfide from lime stabilized biosolids. WEF Odor Conference, Albuquerque, New Mexico (CD-ROM).
29. Oleszkiewicz JA and Mavinic DS (2002). Wastewater biosolids: an overview of processing, treatment, and management. J Environ Eng Sci 1: 75–88. | Article | ChemPort |
30. Papalexandri KP and Dimkou TI (1998). A parametric mixed-integer optimization algorithm for multiobjective engineering problems involving discrete decisions. Ind Eng Chem Res 37: 1866–1882. | Article | ChemPort |
31. ReVelle C and McGarity AE (eds) (1997). Design and Operation of Civil and Environmental Engineering Systems. John Wiley and Sons: New York.
32. Rosenfeld PE (1999). Characterization, quantification, and control of odor emissions from biosolids applications to forest soil. PhD thesis, University of Washington.
33. Ruth JH (1986). Odor thresholds and irritation levels of several chemical substances: a review. Am Ind Hyg Assoc J 47: 142–151.
34. Shih LH and Cheng KJ (2001). Multiobjective transportation planning for waste hauling. J Environ Eng 127: 450–455. | Article | ChemPort |
35. Sostrand P et al (2000). Hazardous peak concentrations of hydrogen sulfide gas related to the sewage purification process. Amer Ind Hyg Assoc J 61: 107–110. | ChemPort |

Appendices

Appendix A: Mathematical details on the Schur decomposition and SOS type 2 variables used

The following lemma is used as part of the approximation scheme to decompose g or (see below for definitions) and involves the determination of the eigenvalues of the Hessian matrix of bilinear functions of the form g or . This result also describes the Schur decomposition of this matrix which is a well-known result in linear algebra (Horn and Johnson, 1985). An advantage of this decomposition is that it involves orthogonal matrices Q which are known to have good numerical stability properties. This lemma will then be used in the main results (Theorem 1). This Appendix and Appendix B are repeated from Gabriel et al (2006) for completeness.

Lemma 1

Let A be a (n+1) (n+1) matrix of the form

where eRⁿis the vector of all ones. Then,

1. The eigenvalues of A are .
2. A=QDQ^T where D=diag(₁..., _n-1, _n, _n+1) and Q is the orthagonal matrix given as
   
   
   
   

Proof

Let eRⁿ be the vector of all ones and partition the vector xRⁿ⁺¹ as x^T=( y^T) with R, yRⁿ. We see that the eigenvalue–eigenvector equations are

For part (i), letting =0 and choosing =0, y0, e^Ty=0 shows that 0 is an eigenvalue of A with geometric multiplicity n-1 so that there are at least n-1 eigenvalues with this value (Horn and Johnson, 1985). Now let , , so that

showing that the other two eigenvalues are .

For part (ii), consider two columns i, j{1, ..., n-1} of the matrix Q with i<j. The inner product of these columns is given by

The inner product with column i and either column n or n+1 is given as

The inner product of columns n and n+1 is given as

The Euclidean norm of column i is i(1/(i+i²))+(i²/(i+i²))=1 and the Euclidean norm for columns n or n+1 is given by (1/2)+n(1/(2n))=1 showing that the matrix is orthagonal. The fact that A=QDQ^T follows easily since the first n-1 columns of Q are eigenvectors associated with a zero eigenvalue in light of the discussion at the start of the proof. Similar reasoning holds for columns n and n+1 which are eigenvectors for the eigenvalues and , respectively. 

This leads to the main result concerning the Hessian matrix of g(x, y)=_i=1ⁿc_ix_iy_i and how to decompose it. In combination with the use of SOS type 2 variable, this theorem justifies the method we propose in this paper.

Theorem 1

Consider the function g: Rⁿ Rⁿ R given by g(x, y)=_i=1ⁿc_ix_iy_i where x, yRⁿ and c_iR₊, i=1..., n. Then, letting z=(x₁ y₁ ... x_n y_n)^T, g(z) is equivalently represented by the separable function (1/2)_i=1ⁿc_iu_i²-(1/2)_i=1ⁿc_iv_i², with the extra linear constraints .

Proof

Letting z=(x₁ y₁ ... x_n y_n)^T, we can write g(z)=(1/2)z^TMz, where M=diag(M₁, ..., M_n) with each

Since M is a real, symmetric matrix, by Schur's decomposition (Horn and Johnson, 1985), there is an orthogonal matrix Q=diag(Q₁, ..., Q_n) such that M=QDQ^T where D=diag(d₁, ..., d_2n). In particular, given the structure of M, from Lemma 1 we have

and (d₁, d₂, ..., d_2n)=(c₁, -c₁, c₂,-c₂, ..., c_n, -c_n). Thus,

with

Remarks

(1) The functions (1/2)c_iu_i², -(1/2)c_iv_i², i=1, ..., n can be approximated with SOS type 2 variables so that g(x, y) can in turn, be approximated with just linear constraints.

(2) The above results can be extended further if one considers N the function : R R^m R defined as . For example, consider the case when n=3. After specifying the order of variables as z=(x₁ y₁ y₂ y₃)^T, the Hessian matrix is

By Lemma 1, ²(z)=QDQ^T, where

Thus,

with

The expression (A.2) can be approximated linearly using SOS type 2 variables.

Appendix B: Approximating the Bilinear Constraints

By Theorem 1 in Appendix A it can be shown that _d=1^|D|(cost _chem)(chem)_d(%flow)_d(B_d)=(1/2)cost _chemd=1^|D|B_d(u_d²-v_d²) with the extra linear constraints , . Next the terms u_d² and v_d² can be piecewise linearized using SOS type 2 variables. Consider the breakpoints b_id^u and b_jd^v for the terms u_d² and v_d², respectively. Note that in practice, it is sometimes computationally advantageous to use the same number of breakpoints for each day as was done for example #2. Then, u_d can be written as _i=1^|I|z_id^ub_id^u and v_d as _j=1^|J|z_jd^vb_jd^v and, u_d², v_d² are replaced respectively, by _i=1^|I|z_id^u(b_id^u)² and _j=1^|J|z_jd^v(b_jd^v)² where _i=1^|I|z_id^u=1, _j=1^|J|z_id^v=1, 0z_id^u1,0z_jd^v1, i=1, ..., |I|, j=1, ..., |J|.

Thus, the bilinear chemical cost term (3) can be approximated by (1/2)cost _chemd=1^|D|B_d(_i=1^|I|z_id^u(b_id^u)²-_j=1^|J|z_jd^v(b_jd^v)², with the additional constraints (XPRESS-MOSEL has SOS type 2 variable types which can be used for a more efficient implementation):

Acknowledgements

We thank the following people for their assistance: Al Razik of MES and Dorian Tolbert formerly of DCWASA.