CN109617048A - A typical scenario selection method for power grid planning based on multi-objective linear programming - Google Patents

Abstract

The invention relates to a method for selecting typical scenarios of power grid planning based on multi-objective linear programming. On the basis of constructing a typical daily evaluation index system, an optimized typical daily selection model and a typical scenario selection method for power grid planning are constructed, including: Evaluate the index system; construct a typical day selection model based on hybrid shaping linear programming. Multi-objective linear programming two-stage fuzzy solution: use the two-stage fuzzy programming solution method to solve the typical day selection model.

Description

A typical scenario selection method for power grid planning based on multi-objective linear programming

technical field

The invention relates to a method for selecting typical scenarios of power grid planning based on multi-objective linear programming.

Background technique

The penetration rate of distributed power in the power grid continues to increase, and the degree of influence on the power grid continues to increase. Therefore, the operation control and planning and design of distribution networks considering various distributed power sources have become a current research focus. However, due to the load point in the power grid Renewable energy sources such as wind and solar have a large amount of data, and the calculation scale required for planning design analysis and scheduling strategy formulation is large. In order to meet the calculation efficiency and the integrity of the data in the calculation, a large amount of data needs to be compressed. selection question

The existing selection of typical daily scenes mainly adopts the clustering method, so more research focuses on the improvement of the clustering algorithm, and combines big data, parallel computing and other technologies to develop efficient algorithms to obtain better results; the other is For a large number of generated scenes, select typical scenes by sampling method; another type is to find typical scenes through heuristic algorithm forward search, backward search and other algorithms, and the improved method is generally optimized by defining new distance indicators and search processes. , to improve the effect of the algorithm, but there is no essential difference between the search and the k-means clustering algorithm. For the selection of typical scenes and typical days, it is necessary to first stipulate how to judge the pros and cons of selection results, and then consider how to select them. However, the existing literature does not solve the problem of evaluation indicators, and almost all use distance indicators to evaluate selection results.

In recent years, some research results have been achieved on the selection of typical scenes and typical days, but the existing methods still have certain defects and deficiencies:

(1) None of the existing literatures have solved the problem of evaluation indicators, and only the distance indicators are used to evaluate the selection results, which lacks a comprehensive basis.

(2) Most of the existing selection methods for typical daily scenes use the clustering method, so more researches focus on the improvement of the clustering algorithm.

SUMMARY OF THE INVENTION

In view of the above problems, the purpose of the present invention is to overcome the deficiencies of a single typical day selection method, and on the basis of constructing a typical day evaluation index system, construct an optimal typical day selection model and selection algorithm, and the technical scheme is as follows:

A method for selecting typical scenarios of power grid planning based on multi-objective linear programming. On the basis of constructing a typical daily evaluation index system, an optimal typical daily selection model and a method for selecting typical scenarios of power grid planning are constructed, including:

S1) Typical daily evaluation index system

1) Statistical indicators

The annual total load power deviation ΔC represents the relative error between the total load power ∑ω _d ·C _d and the total load power C _year of the original data after weighted calculation on a typical day:

In the above formula, ω _d represents the weight coefficient of the typical day d, C _d represents the total electricity load of the typical day, C _year represents the total load electricity of the year, and D represents the set of all typical days;

The annual load power distribution deviation ΔP represents the total load power calculated by weighting for a typical day in each period and the total historical load at this moment The mean relative error of :

In the above formula, D ₀ represents the set of all historical dates in the original data, represents the original load power value of date d at time t, Represents the load power value at the t-th time on a typical day;

The annual resource total deviation ΔS represents the relative error between the total resource amount ∑ω _d ·S _d and the total resource amount S _year in the original data after weighted calculation on a typical day; where S _d represents the total resource amount on a typical day d, S _year Represents the total amount of resources for the year;

The annual resource distribution deviation ΔW represents the total resource amount after weighted calculation for typical days in each period and the total amount of historical resources at this moment The mean relative error of ; where, represents the original resource value of date d at time t, Represents the resource value at the t-th time on a typical day;

2) Timing indicators

The data density around a typical day is represented by the number of data points within the truncation distance:

I _S = {1,2,...,card(D ₀ )}

In the above formula, d _ij represents the distance between the i-th and j-th typical daily data vectors respectively, the present invention adopts the Euclidean distance, d _c represents the truncation distance, and _IS represents the index set;

The typical daily radiation radius is defined by distance. If the typical daily i is the data point with the largest data density in the world, the radiation radius is the distance between the point and the global farthest point, otherwise it is determined as the distance between the nearest adjacent data point with greater data density. spacing:

In the above formula, Represents the index set of the ith typical day, which is composed of the labels of individuals with larger surrounding data density;

The peak load deviation ΔL represents the typical mid-day maximum load value at the same time The maximum load value at the corresponding moment in the historical data relative error;

The peak resource deviation ΔS represents the maximum resource value in a typical day at the same time The maximum resource value at the corresponding moment in the historical data relative error;

Maximum deviation of power change rate by time period Reflects the maximum load variation power during a certain period of time in a typical day and the largest change in historical data relative error;

Maximum deviation of resource change rate by time period Reflects the maximum resource change value for a typical day in a certain time period and the largest change in historical data relative error;

Coverage of power change rate by time period Reflects the maximum load fluctuation power for a certain period of time in a typical day Relative position in historical data change value:

Time-by-period resource change rate coverage Reflects the maximum resource change value for a certain period of time in a typical day. The relative position in the historical data change value, the composition form is the same as the time period power change rate coverage

S2), build a typical day selection model based on hybrid shaping linear programming

The optimization objective z ₁ represents the selected number of typical days, z ₂ represents the error of the total load demand and total resource amount calculated by weighting on the typical day, z ₃ represents the data density around the total typical day, and z ₄ represents the total typical day. Radiation Radius:

In the above formula, u _i represents the binary variable selected for a typical day, u _i is 1, it means that the ith day is a typical day, n represents the total number of days, ρ=[ρ ₁ ,ρ ₂ ,...,ρ _n ]δ=[δ ₁ ,δ ₂ ,…,δ _n ], the column in matrix A represents the load data and resource data of one day, and b represents the total amount of load and resource at the same time:

The optimization variables include the weight variable _wi and the binary variable _ui :

In the constraints (1) the weight of typical days is constrained by binary variables, if it is not a typical day, the weight is reset to zero; (2) the sum of the weights of all typical days is the total number of days N in the historical data; (3) the number of days in each period is represented. Load or resource deviation should keep the total deviation within a certain range, and α is the proportional coefficient; (4) Set the lower limit of typical days, and extreme constraints can be set by constraining the formulation time; (5) Indicates that the typical daily weight is not Negative real number; (6) indicates that the variable _ui is a binary variable;

S3) Two-stage fuzzy solution of multi-objective linear programming

The typical day selection model is solved by a two-stage fuzzy programming method.

The present invention has the following advantages compared with the prior art due to the adoption of the above technical solutions:

(1) The present invention firstly constructs an evaluation index system for evaluating the quality of the selection results of typical days, and then constructs a hybrid shaping linear programming model according to the evaluation index, and optimizes the allocation weights, so that in the case of selecting the minimum number of days, resources and load The total deviation and distribution error are the smallest, and the typical day has a good representativeness.

(2) The optimization method proposed by the present invention can solve the optimization result of the weight, so the error in the total amount and distribution deviation is smaller than that of the traditional clustering algorithm, and the linear programming model can flexibly set the extreme value deviation, fluctuation deviation and Various extreme constraints such as typical number of days make the selection result in line with the expected effect and usage.

Description of drawings

Fig. 1 is a typical daily selection comprehensive evaluation index structure diagram of the present invention.

Figure 2 is a month of resource and load data for an embodiment.

FIG. 3 is a decision diagram based on density and radiation radius of an embodiment.

FIG. 4 is a diagram of a situation of raw data corresponding to a typical data center of the embodiment.

Fig. 5 is the selection result decision position diagram of embodiment MILP

Fig. 6 is the typical day selection result comparison chart of the embodiment

FIG. 7 is a graph of peak load deviation by time period of an embodiment

FIG. 8 is a statistic diagram of peak resource deviation by time period according to an embodiment

Description of the labels in the figure: 1 black: k-means method; 2 red line: MILP method; 3 blue line: MILP2 method

Table 1 shows the extreme values of the respective objectives in the single-objective planning of the embodiment of the present invention.

Table 2 shows the dates when each time period in the embodiment of the present invention satisfies the deviation threshold.

Table 3 is a comparison of the solution results with extreme constraints and without extreme constraints in the embodiment of the present invention.

Table 4 is a comparison of various indicators of the MILP and k-means method with extreme constraints and without extreme constraints in the embodiment of the present invention.

Detailed ways

There is a large amount of information in historical data, including total resource and load information, distribution information, time series change information, and extreme scene information. A typical daily comprehensive evaluation index model including statistical indicators and time series indicators is constructed to comprehensively evaluate from various aspects. Selection effect of a typical day. A multi-objective hybrid shaping linear programming model is constructed by combining statistical indicators and time-series indicators, in order to reduce the number of typical days and optimize the indicators as the goal of building the model for the selection of typical days.

On the basis of constructing a typical day evaluation index system, the present invention constructs an optimal typical day selection model and a selection algorithm. The present invention will be described in detail below with reference to the accompanying drawings and embodiments.

1. Implementation method

1) Comprehensive evaluation index system

A large amount of information in historical data includes total resource and load information, distribution information, timing change information, and extreme scene information. Therefore, in order to comprehensively evaluate the selection effect of typical days from various aspects, the present invention constructs a comprehensive evaluation index model for typical day selection including statistical indicators and time series indicators. The model structure diagram is shown in FIG. 1 .

The establishment of statistical indicators mainly considers that because the total amount and distribution of loads and resources involve the calculation of investment benefits of distributed power and distribution networks, the statistical indicators of historical data and typical daily data should be kept within a certain error range.

The annual total load power deviation ΔC represents the relative error between the total load power ∑ω _d ·C _d and the total load power C _year of the original data after weighted calculation on a typical day:

In the above formula, ω _d represents the weight coefficient of the typical day d, C _d represents the total electricity load of the typical day, C _year represents the total load electricity of the year, and D represents the set of all typical days.

The annual load power distribution deviation ΔP represents the total load power calculated by weighting for a typical day in each period and the total historical load at this moment The mean relative error of :

In the above formula, D ₀ represents the set of all historical dates in the original data, represents the original load power value of date d at time t, Indicates the load power value at the t-th time on a typical day.

The annual resource total deviation ΔS represents the relative error between the total resource amount Σω _d ·S _d after weighted calculation on a typical day and the total resource amount S _year in the original data. Among them, S _d represents the total resources of a typical day d, and S _year represents the total resources of the whole year.

The annual resource distribution deviation ΔW represents the total resource amount after weighted calculation for typical days in each period and the total amount of historical resources at this moment The mean relative error of . in, represents the original resource value of date d at time t, Represents the resource value at the t-th time on a typical day.

Secondly, for the typical day selection of a large amount of time series data, it is necessary to select a scenario with a high typical degree and often occurs in actual operation. Therefore, the present invention constructs an index reflecting the typical day typical degree. The typical day surrounding data density and typical daily radiation radius. In addition, the establishment of time sequence indicators should also consider the existence of holidays, bad weather, etc. in actual operation. Therefore, the present invention constructs the deviation between peak resources and peak load to reflect the description effect of typical days on extreme scenarios. Due to the time-series fluctuation of resources and loads, which affects the operating frequency of the distribution network, the ramping of units, and the maximum charging and discharging power of the energy storage device, the volatility reflected in a typical day and the extreme volatility in the original data should be maintained. within a certain error range.

The data density around a typical day is represented by the number of data points within the truncation distance:

I _S = {1,2,...,card(D ₀ )}

In the above formula, d _ij represents the distance between the i-th and j-th typical day data vectors respectively, the present invention adopts the Euclidean distance, d _c represents the cut-off distance, and _IS represents the index set.

The typical daily radiation radius is defined by distance. If the typical daily i is the data point with the largest data density in the world, the radiation radius is the distance between the point and the global farthest point, otherwise it is determined as the distance between the nearest adjacent data point with greater data density. spacing:

In the above formula, The set of indicators representing the ith typical day, consisting of the labels of other individuals whose surrounding data density is greater than itself.

The peak load deviation ΔL represents the typical mid-day maximum load value at the same time The maximum load value at the corresponding moment in the historical data relative error.

The peak resource deviation ΔS represents the maximum resource value in a typical day at the same time The maximum resource value at the corresponding moment in the historical data relative error.

Maximum deviation of power change rate by time period Reflects the maximum load variation power during a certain period of time in a typical day and the largest change in historical data relative error.

Maximum deviation of resource change rate by time period Reflects the maximum resource change value for a typical day in a certain time period and the largest change in historical data relative error.

Coverage of power change rate by time period Reflects the maximum load fluctuation power for a certain period of time in a typical day Relative position in historical data change value:

Time-by-period resource change rate coverage Reflects the maximum resource change value for a certain period of time in a typical day. The relative position in the historical data change value, the composition form is the same as the time period power change rate coverage

2) Typical day selection model

The main purpose of the typical day selection is to select as few typical days as possible to replace a large amount of original data. Therefore, in the case of reducing the number of typical days as much as possible, making the optimal indicators become the basis for the construction of the typical day selection model. The above statistical indicators and time series indicators construct a multi-objective mixed shaping linear programming model.

The optimization objective z ₁ represents the selected number of typical days, z ₂ represents the error of the total load demand and total resource amount calculated by weighting on the typical day, z ₃ represents the data density around the total typical day, and z ₄ represents the total typical day. Radiation Radius:

In the above formula, u _i represents the binary variable selected for a typical day, u _i is 1, it means that the ith day is a typical day, n represents the total number of days, ρ=[ρ ₁ ,ρ ₂ ,...,ρ _n ]δ=[δ ₁ ,δ ₂ ,…,δ _n ], each column of the A matrix represents the load data and resource data of each day, and b represents the total amount of load and resources at the same time:

The optimization variables include the weight variable _wi and the binary variable _ui :

In the constraints (1) the weight of typical days is constrained by binary variables, if it is not a typical day, the weight is reset to zero; (2) the sum of the weights of all typical days is the total number of days N in the historical data; (3) the number of days in each period is represented. Load or resource deviation should keep the total deviation within a certain range, and α is the proportional coefficient; (4) Set the lower limit of typical days, and extreme constraints can be set by constraining the formulation time; (5) Indicates that the typical daily weight is not Negative real number; (6) indicates that the variable _ui is a binary variable

3) Multi-objective linear programming solution

The above multi-objective linear programming model is solved by a two-stage fuzzy programming method.

2. Case study

1) Validation of typical indicators

The case data uses the one-month wind speed, light intensity and local electricity load measured in a certain place in Xining. The data after the per-unitization of resources and loads are shown in Figure 2.

In order to verify the validity of the proposed indicators of typical daily surrounding data density and typical daily radiation radius, the following schematic diagram is constructed with the data density as the abscissa and the radiation radius as the ordinate, as shown in Figure 3. The point set near the upper right corner of the figure reflects the The density of the output data is large and the radiation radius is large, which is the center point of a large amount of data, and the typicality is strong. Therefore, p1 and p2 in the upper right corner are selected as representatives; the point in the upper left corner indicates that the surrounding density value of the data is small, and the radiation radius is small. Large, representing outliers in extreme cases, so select the upper left corner p3 and p4 for comparison. The original data graph for plotting p1-p4 is shown in Figure 4.

It can be seen from Figure 5 that the data density around a typical day and the typical daily radiation radius can better reflect the typical pattern. The wind speed, light intensity, and load demand in the typical data p1 and p2 are all around the average value of the overall data, while the wind speed in p3 and p4, which represent extreme scenarios, is high and the light resources are insufficient, indicating that the weather conditions on the day are not good. It is cloudy or rainy. The selected typical days cover a variety of states, so the typicality indicator can better reflect the typicality of the data.

2) Selection of typical days without extreme constraints

Table 1 shows that the maximum reference real value of the objective function is set as the expected value, and the maximum value of the typical days of Obj1 is set as the boundary value.

Table 3 shows the two-stage solution results of the fuzzy programming method in the present invention under two conditions of the extreme constraint of the minimum value of the air resource and the extreme constraint of the minimum value of the air resource not added. If the two-stage solution results are the same, the results are valid. The most typical data on the 12th day is given the largest weight value to ensure the rationality and representativeness of the selection results, while the less typical data on the 73rd and 90th days are given a smaller value. Weight, the selection result includes both typical data and atypical data, and the practicability and rationality of the typical day selection result is ensured through reasonable weight distribution. Typical daily results selected in the example MILP are shown in FIG. 5 .

Table 4 shows the comparison of various indicators between the typical daily method and the k-means method selected by the multi-objective mixed shaping linear programming (MILP) method in the present invention. It can be seen that in terms of load and total resources, since the weight coefficients of the clustering method are all positive integers, and the method proposed by the present invention optimizes the real value of the weight, the error in the total statistics is smaller, and the accuracy can generally be improved by 10 times. about. By setting the error constraint of each time period, it can be seen that the distribution error of resource and load is also smaller than the result of k-means method. There are no mandatory constraints on the indicators such as resource and load fluctuation coverage and extreme maximum deviation. Due to the constraints of the total deviation by time period, there are some extreme scenarios in the selection results. It can be seen that the deviation of the MILP method other than load is relatively In addition, in terms of resource and load coverage, except for the fluctuation rate of light coverage, which is basically the same as that of the k-means method, it is superior in the peak deviation of resource and load fluctuation, indicating that the extreme value of the selection results on a typical day is more than the extreme value of the original data. Closer, it can better reflect the extreme situation that actually exists.

3) Add extreme constraints to select typical days

Under the condition of no extreme constraints, the deviation of the minimum value of wind speed is too large, and extreme constraints of the minimum value of wind speed can be added. Since there are sufficient light resources during the period from 11:00 to 13:00, the size of wind resources mainly depends on the operation of the power grid and energy storage. Therefore, it is stipulated that the relative deviation of the minimum wind speed during this period is less than 1. Table 2 gives the date numbers that satisfy the set thresholds corresponding to the three times respectively.

Table 3 shows the solution results of the two-stage planning after adding the wind speed minimum extreme constraint. After considering the extreme index, the most representative day 12 still occupies the main weight, and the extreme of the 58th day satisfies the constraint of the minimum wind speed. The planning result includes this day and assigns a smaller weight. In Table 4, the decrease of the typicality index without extreme constraint conditions further reduces the total deviation and distribution deviation, all of which are better than the k-means algorithm. The coverage of lighting resources is further improved, and the fluctuation peak deviation is further reduced.

Figure 6 shows the original data of the k-means algorithm and the MILP2 selection result after adding extreme constraints. It can be seen that the extreme constraints of the present invention make the distribution of wind resources wider, and the scenarios of the maximum wind resources and the minimum wind resources are the same as those in the original data. The extreme values are closer.

Figure 7 compares the deviation of the peak load by time period of the three methods. MILP2 represents the case of adding wind speed constraints. It can be seen that the peak load deviation after adding wind speed constraints is reduced, but in general, the three situations are basically the same, and all of them can make the load The peak error is kept within a small range.

FIG. 8 compares the deviation of peak resources by time period in three cases. The deviation between the maximum and minimum values of visible light resources can be obtained through MILP2 of the method of the present invention, which is smaller than that of k-means. After adding the wind speed minimum constraint, it can be seen from Figure 8 (4) that the result of the selection greatly reduces the deviation of the minimum wind resource value. The deviation values at the three noon time are all within 1, and the errors in other aspects are basically maintained compared with those without constraints. constant.

Table 1

Target max min Obj1 Typical days of the day 15 3 Obj2 Total Resource/Load Deviation Rate 0 0.1539 Obj3 total density 2022.6740 256.4277 Obj4 total radiation radius 171.4572 23.9521

Table 2

table 3

Table 4

Various indicators k-means MILP MILP2 Load capacity deviation 0.017084 0.001108 0.00047 Light resource deviation 0.067755 0.002045 0.000186 Wind resource deviation 0.074744 0.006427 0.002769 Load distribution error 0.011262 0.017118 0.00963 Light resource distribution error 0.067755 0.034254 0.038195 Wind resource distribution error 0.080398 0.032887 0.039933 load fluctuation coverage 0.988406 0.999034 0.999034 Load fluctuation peak deviation 0.207638 0.039367 0.039367 Light Fluctuation Coverage 0.997585 0.987923 0.993237 Light fluctuation peak deviation 0.04228 0.190987 0.139173 Wind Speed Fluctuation Coverage 0.987923 0.994203 0.994203 Wind speed fluctuation peak deviation 0.588168 0.509592 0.509592

Claims (1)

1. A power grid planning typical scene selection method based on multi-objective linear programming is used for constructing an optimized typical day selection model and a power grid planning typical scene selection method on the basis of constructing a typical day evaluation index system, and comprises the following steps:

s1) typical daily evaluation index system

1) Statistical index

The annual total load electric quantity deviation delta C represents the total load electric quantity sigma omega after the typical day is calculated through weighting_d·C_dTotal load capacity C with original data_yearRelative error of (2):

in the above formula, ω_dWeight coefficient, C, representing typical day d_dTotal electrical load capacity of the whole day, C, representing typical day d_yearRepresenting the total annual load capacity, and D representing the set of all typical days;

the annual load power distribution deviation Δ P represents the total load capacity calculated by weighting for each period of a typical dayAnd the total amount of the historical load at the momentAverage value of relative error of (a):

in the above formula, D₀Representing the set of all historical dates in the raw data,indicating the original load power value at time t on date d,representing the load power value at the tth moment of a typical day d;

the annual resource total deviation Delta S represents the total resource amount Sigma omega after the typical day is calculated by weighting_d·S_dAnd the total amount S of resources in the original data_yearRelative error of (2); wherein S_dRepresents the total amount of resources, S, for a typical day d_yearRepresenting the total annual resource amount;

the annual resource distribution deviation aw represents the total resource amount calculated by weighting for each period of the typical dayAnd the total amount of historical resources at the momentAverage value of relative error of; wherein,the raw asset value representing the date d at time t,a resource value representing the typical day, d, time t;

2) timing indicator

Typical day-around data density is represented by the number of data points within a cutoff distance:

I_S＝{1,2,…,card(D₀)}

in the above formula, d_ijRespectively representing the distance between the ith and jth typical day data vectors, the Euclidean distance, d, is adopted in the invention_cDenotes the truncation distance, I_SRepresenting a set of metrics;

the typical daily radiance radius is defined by using distance, if the typical day i is the global maximum data density data point, the radiance radius is the distance between the point and the global farthest point, otherwise, the radiance radius is the distance between the point and the adjacent closest data point with greater data density:

in the above formula, the first and second carbon atoms are,the index set representing the ith typical day is composed of labels of individuals with higher surrounding data density;

the peak load deviation Δ L represents the maximum load value in a typical day at the same timeMaximum load value at time corresponding to historical dataRelative error of (2);

the peak resource deviation Delta S represents the maximum resource value in a typical day at the same timeMaximum resource value corresponding to time in historical dataRelative error of (2);

maximum deviation of time-interval power change rateReflecting the maximum load change power in a certain period of time in a typical dayAnd the maximum change value in the historical dataRelative error of (2);

maximum deviation of time-interval resource change rateTypical day of the reactionMaximum resource change value of a certain period of timeAnd the maximum change value in the historical dataRelative error of (2);

power change rate coverage over time periodsReflecting the maximum load variation power for a certain period of time on a typical dayRelative position in historical data variance:

time-phased resource change rate coverageReflecting the maximum resource variation value for a certain period of time in a typical day.Relative position in historical data variation value, form same time period power change rate coverage

S2), constructing a typical day selection model based on mixed shaping linear programming

Optimization objective z₁Represents typical days of selection, z₂Error, z, representing the total load demand and total resource volume after a typical day has been calculated by weighting₃Representing the total typical data density around the day, z₄Represents the total typical daily radiation radius:

in the above formula, u_iBinary variable, u, representing typical day picks_i1 denotes that the i-th day is a typical day, n denotes the total number of days, and ρ ═ ρ₁,ρ₂,…,ρ_n]δ＝[δ₁,δ₂,…,δ_n]The columns in the matrix A represent load data and resource data of one day, and b represents the total amount of load and resource at the same time:

the optimization variables include a weight variable w_iAnd binary variable u_i：

The constraint conditions comprise (1) constraint of typical day weights through binary variables, and zero setting of the weights if the typical day weights are not, 2) representation of sum of all typical day weights as total days N in historical data, (3) representation of load or resource deviation of each time interval to enable the total deviation to be controlled within a certain range, α is a proportionality coefficient, (4) setting of lower limit of typical day days, extreme constraint can be set through constraint making time, (5) representation of non-negative real numbers of typical day weights, and (6) representation of variable u_iIs a binary variable;

s3) Multi-object Linear programming two-stage fuzzy solution

And solving the typical day selection model by adopting a two-stage fuzzy programming solution.