www.gusucode.com > optim 案例源码 matlab代码程序 > optim/airpollution.m

    %% Analyzing the Effect of Uncertainty Using Semi-Infinite Programming
% 
% This example shows how to use semi-infinite programming to investigate
% the effect of uncertainty in the model parameters of an optimization
% problem. We will formulate and solve an optimization problem using the
% function |fseminf|, a semi-infinite programming solver in Optimization Toolbox(TM).
% 
% The problem illustrated in this example involves the control of air
% pollution. Specifically, a set of chimney stacks are to be built in a
% given geographic area. As the height of each chimney stack increases, the
% ground level concentration of pollutants from the stack decreases.
% However, the construction cost of each chimney stack increases with
% height. We will solve a problem to minimize the cumulative height of the
% chimney stacks, hence construction cost, subject to ground level
% pollution concentration not exceeding a legislated limit. This problem is
% outlined in the following reference:
% 
% Air pollution control with semi-infinite programming, A.I.F. Vaz and E.C.
% Ferreira, XXVIII Congreso Nacional de Estadistica e Investigacion
% Operativa, October 2004
%
% In this example we will first solve the problem published in the above
% article as the _Minimal Stack Height_ problem. The models in this problem
% are dependent on several parameters, two of which are wind speed and
% direction. All model parameters are assumed to be known exactly in the
% first solution of the problem.
%
% We then extend the original problem by allowing the wind speed and
% direction parameters to vary within given ranges. This will allow us to
% analyze the effects of uncertainty in these parameters on the optimal
% solution to this problem.

%   Copyright 2008-2012 The MathWorks, Inc.

%% Minimal Stack Height Problem
% Consider a 20km-by-20km region, R, in which ten chimney stacks are to be
% placed. These chimney stacks release several pollutants into the
% atmosphere, one of which is sulfur dioxide. The x, y locations of the
% stacks are fixed, but the height of the stacks can vary. 
%
% Constructors of the stacks would like to minimize the total height of the
% stacks, thus minimizing construction costs. However, this is balanced by
% the conflicting requirement that the concentration of sulfur dioxide at
% any point on the ground in the region R must not exceed the legislated
% maximum.

%%
% First, let's plot the chimney stacks at their initial height. Note that
% we have zoomed in on a 4km-by-4km subregion of R which contains the
% chimney stacks.
h0 = [210;210;180;180;150;150;120;120;90;90];
plotChimneyStacks(h0, 'Chimney Stack Initial Height');

%%
% There are two environment related parameters in this problem, the wind
% speed and direction. Later in this example we will allow these parameters to
% vary, but for the first problem we will set these parameters to typical
% values.

% Wind direction in radians
theta0 = 3.996;   
% Wind speed in m/s
U0 = 5.64;    

%%
% Now let's plot the ground level concentration of sulfur dioxide (SO2)
% over the entire region R (remember that the plot of chimney stacks was
% over a smaller region). The SO2 concentration has been calculated with
% the chimney stacks set to their initial heights. 
%
% We can see that the concentration of SO2 varies over the region of
% interest. There are two features of the Sulfur Dioxide graph of note:
% 
% * SO2 concentration rises in the top left hand corner of the (x,y) plane
% * SO2 concentration is approximately zero throughout most of the region
% 
% In very simple terms, the first feature is due to the prevailing wind,
% which is blowing SO2 toward the top left hand corner of the (x,y) plane
% in this example. The second factor is due to SO2 being transported to the
% ground via diffusion. This is a slower process compared to the prevailing
% wind and thus SO2 only reaches ground level in the top left hand corner
% of the region of interest.
%
% For a more detailed discussion of atmospheric dispersion from chimney
% stacks, consult the reference cited in the introduction.
%
% The pink plane indicates a SO2 concentration of $$ 0.000125 gm^{-3} $$.
% This is the legislated maximum for which the Sulfur Dioxide concentration
% must not exceed in the region R. It can be clearly seen from the graph
% that the SO2 concentration exceeds the maximum for the initial chimney
% stack height. 
%
% Examine the MATLAB file |concSulfurDioxide| to see how the sulfur dioxide
% concentration is calculated.

plotSulfurDioxide(h0, theta0, U0, ...
    'Sulfur Dioxide Concentration at Initial Stack Height');

%% How |fseminf| Works
%
% Before we solve the minimal stack height problem, we will outline how
% |fseminf| solves a semi-infinite problem. A general semi-infinite
% programming problem can be stated as:
%
% $$\min f(x)$$
%
% such that 
% 
% $$Ax <= b $$ (Linear inequality constraints)
%
% $$Aeq*x = beq $$ (Linear equality constraints)
%
% $$c(x) <= 0 $$ (Nonlinear Inequality Constraints)
%
% $$ceq(x) = 0 $$ (Nonlinear Equality Constraints)
%
% $$l <= x <= u $$ (Bounds)
%
% and 
%
% $$K_{j}(x, w) <= 0 $$, where $$w \in I_{j} $$ for $$j = 1,..., n_{inf} $$
% (Nonlinear semi-infinite constraints)
%
% This algorithm allows you to specify constraints for a nonlinear
% optimization problem that must be satisfied over intervals of an
% auxiliary variable, $$ w $$. Note that for |fseminf|, this variable is
% restricted to be either 1 or 2 dimensional for each semi-infinite
% constraint.
%
% The function |fseminf| solves the general semi-infinite problem by
% starting from an initial value, $$ x_{0} $$, and using an iterative
% procedure to obtain an optimum solution, $$ x_{opt} $$. 
% 
% The key component of the algorithm is the handling of the "semi-infinite"
% constraints, $$ K_{j} $$. At $$ x_{opt} $$ it is required that the $$
% K_{j} $$ must be feasible at every value of $$ w $$ in the interval $$
% I_{j} $$. This constraint can be simplified by considering all the local 
% maxima of $$ K_{j} $$ with respect to $$ w $$ in the interval $$ I_{j}
% $$. The original constraint is equivalent to requiring that the value of
% $$ K_{j} $$ at each of the above local maxima is feasible.
%
% |fseminf| calculates an approximation to all the local maximum values of
% each semi-infinite constraint, $$ K_{j} $$. To do this, |fseminf| first
% calculates each semi-infinite constraint over a mesh of $$ w $$ values. A
% simple differencing scheme is then used to calculate all the local
% maximum values of $$ K_{j} $$ from the evaluated semi-infinite
% constraint.
%
% As we will see later, you create this mesh in your constraint function.
% The spacing you should use for each $$ w $$ coordinate of the mesh is
% supplied to your constraint function by |fseminf|.
%
% At each iteration of the algorithm, the following steps are performed:
% 
% # Evaluate $$ K_{j} $$ over a mesh of $$ w $$-values using the current
% mesh spacing for each $$ w $$-coordinate.
% # Calculate an approximation to all the local maximum values of $$ K_{j}
% $$ using the evaluation of $$ K_{j} $$ from step 1.
% # Replace each $$ K_{j} $$ in the general semi-infinite problem with the
% set of local maximum values found in steps 1-2. The problem now has a
% finite number of nonlinear constraints. |fseminf| uses the SQP algorithm
% used by |fmincon| to take one iteration step of the modified problem.
% # Check if any of the SQP algorithm's stopping criteria are met at the
% new point $$ x $$. If any criteria are met the algorithm terminates; if
% not, |fseminf| continues to step 5. For example, if the first order
% optimality value for the problem defined in step 3 is less than the
% specified tolerance then |fseminf| will terminate.
% # Update the mesh spacing used in the evaluation of the semi-infinite
% constraints in step 1.

%% Writing the Nonlinear Constraint Function
%
% Before we can call |fseminf| to solve the problem, we need to write a
% function to evaluate the nonlinear constraints in this problem. The
% constraint to be implemented is that the ground level Sulfur Dioxide
% concentration must not exceed $$ 0.000125 gm^{-3} $$ at every point in
% region R.
%
% This is a semi-infinite constraint, and the implementation of the
% constraint function is explained in this section. For the minimal stack
% height problem we have implemented the constraint in the MATLAB file
% |airPollutionCon|.

type airPollutionCon.m

%%
% This function illustrates the general structure of a constraint function
% for a semi-infinite programming problem. In particular, a constraint
% function for |fseminf| can be broken up into three parts:
% 
% 1. Define the initial mesh size for the constraint evaluation
%
% Recall that |fseminf| evaluates the "semi-infinite" constraints over a
% mesh as part of the overall calculation of these constraints. When your
% constraint function is called by |fseminf|, the mesh spacing you should
% use is supplied to your function. |Fseminf| will initially call your
% constraint function with the mesh spacing, |s|, set to NaN. This
% allows you to initialize the mesh size for the constraint evaluation.
% Here, we have one "infinite" constraint in two "infinite" variables. This
% means we need to initialize the mesh size to a 1-by-2 matrix, in this
% case, |s = [1000 4000]|.
%
% 2. Define the mesh that will be used for the constraint evaluation
%
% A mesh that will be used for the constraint evaluation needs to be
% created. The three lines of code following the comment "Define the grid
% that the "infinite" constraints will be evaluated over" in
% |airPollutionCon| can be modified for most 2-d semi-infinite programming
% problems.
%
% 3. Calculate the constraints over the mesh
%
% Once the mesh has been defined, the constraints can be calculated over
% it. These constraints are then returned to |fseminf| from the above
% constraint function.
%
% Note that in this problem, we have also rescaled the constraints so that
% they vary on a scale which is closer to that of the objective function.
% This helps |fseminf| to avoid scaling issues associated with objectives
% and constraints which vary on disparate scales.

%% Solve the Optimization Problem
%
% We can now call |fseminf| to solve the problem. The chimney stacks must
% all be at least 10m tall and we use the initial stack height specified
% earlier. Note that the third input argument to |fseminf| below (1)
% indicates that there is only one semi-infinite constraint.

lb = 10*ones(size(h0));
[hsopt, sumh, exitflag] = fseminf(@(h)sum(h), h0, 1, ...
    @(h,s) airPollutionCon(h,s,theta0,U0), [], [], [], [], lb);
fprintf('\nMinimum computed cumulative height of chimney stacks : %7.2f m\n', sumh);

%% 
% The minimum cumulative height computed by |fseminf| is considerably
% higher than the initial total height of the chimney stacks. We will see
% how the minimum cumulative height changes when parameter uncertainty is
% added to the problem later in the example. For now, let's plot the chimney
% stacks at their optimal height.
%
% Examine the MATLAB file |plotChimneyStacks| to see how the plot was generated.

plotChimneyStacks(hsopt, 'Chimney Stack Optimal Height');

%% Check the Optimization Results
%
% Recall that |fseminf| determines that the semi-infinite constraint is
% satisfied everywhere by ensuring that discretized maxima of the
% constraint are below the specified bound. We can verify that the
% semi-infinite constraint is satisfied everywhere by plotting the ground
% level sulfur dioxide concentration for the optimal stack height.
%
% Note that the sulfur dioxide concentration takes its maximum possible
% value in the upper left corner of the (x, y) plane, i.e. at x = -20000m,
% y = 20000m. This point is marked by the blue dot in the figure below and
% verified by calculating the sulfur dioxide concentration at this point.
%
% Examine the MATLAB file |plotSulfurDioxide| to see how the plots was
% generated.

titleStr = 'Optimal Sulfur Dioxide Concentration and its maximum (blue)';
xMaxSD = [-20000 20000];
plotSulfurDioxide(hsopt, theta0, U0, titleStr, xMaxSD);
SO2Max = concSulfurDioxide(-20000, 20000, hsopt, theta0, U0);
fprintf('Sulfur Dioxide Concentration at x = -20000m, y = 20000m : %e g/m^3\n', SO2Max);  

%% Considering Uncertainty in the Environmental Factors
% The sulfur dioxide concentration depends on several environmental factors
% which were held at fixed values in the above problem. Two of the
% environmental factors are wind speed and wind direction. See the
% reference cited in the introduction for a more detailed discussion of all
% the problem parameters.
%
% We can investigate the change in behavior for the system with respect to
% the wind speed and direction. In this section of the example, we want to
% make sure that the sulfur dioxide limits are satisfied even if the wind
% direction changes from 3.82 rad to 4.18 rad and mean wind speed varies
% between 5 and 6.2 m/s.

%%
% We need to implement a semi-infinite constraint to ensure that the sulfur
% dioxide concentration does not exceed the limit in region R. This
% constraint is required to be feasible for all pairs of wind speed and
% direction.
% 
% Such a constraint will have four "infinite" variables (wind speed and
% direction and the x-y coordinates of the ground). However, any
% semi-infinite constraint supplied to |fseminf| can have no more than two
% "infinite" variables.
%
% To implement this constraint in a suitable form for |fseminf|, we recall
% the SO2 concentration at the optimum stack height in the previous
% problem. In particular, the SO2 concentration takes its maximum possible
% value at x = -20000m, y = 20000m. To reduce the number of "infinite"
% variables, we will assume that the SO2 concentration will also take its
% maximum value at this point when uncertainty is present. We then require
% that SO2 concentration at this point is below $$ 0.000125 gm^{-3} $$ for
% all pairs of wind speed and direction.
%
% This means that the "infinite" variables for this problem are wind speed
% and direction. To see how this constraint has been implemented, inspect
% the MATLAB file |uncertainAirPollutionCon|.

type uncertainAirPollutionCon.m

%% 
% This constraint function can be divided into same three sections as
% before:
%
% 1. Define the initial mesh size for the constraint evaluation
%
% The code following the comment "Initial sampling interval" initializes
% the mesh size.
%
% 2. Define the mesh that will be used for the constraint evaluation
%
% The next section of code creates a mesh (now in wind speed and direction)
% using a similar construction to that used in the initial problem.
%
% 3. Calculate the constraints over the mesh
% 
% The remainder of the code calculates the SO2 concentration at each point
% of the wind speed/direction mesh. These constraints are then returned to
% |fseminf| from the above constraint function.

%%
% We can now call |fseminf| to solve the stack height problem considering
% uncertainty in the environmental factors. 

[hsopt2, sumh2, exitflag2] = fseminf(@(h)sum(h), h0, 1, ...
    @uncertainAirPollutionCon, [], [], [], [], lb);
fprintf('\nMinimal computed cumulative height of chimney stacks with uncertainty: %7.2f m\n', sumh2);

%%
% We can now look at the difference between the minimum computed cumulative
% stack height for the problem with and without parameter uncertainty. You
% should be able to see that the minimum cumulative height increases when
% uncertainty is added to the problem. This expected increase in height
% allows the SO2 concentration to remain below the legislated maximum for
% all wind speed/direction pairs in the specified range.
%
% We can check that the sulfur dioxide concentration does not exceed the
% limit over the region of interest via inspection of a sulfur dioxide
% plot. For a given (x, y) point, we plot the maximum SO2
% concentration for the wind speed and direction in the stated ranges. Note
% that we have zoomed in on the upper left corner of the X-Y plane.
titleStr = 'Optimal Sulfur Dioxide Concentration under Uncertainty';
thetaRange = 3.82:0.02:4.18;
URange = 5:0.2:6.2;
XRange = [-20000,-15000];
YRange = [15000,20000];
plotSulfurDioxideUncertain(hsopt2, thetaRange, URange, XRange, YRange, titleStr);

%% 
% We finally plot the chimney stacks at their optimal height when there is 
% uncertainty in the problem definition.

plotChimneyStacks(hsopt2, 'Chimney Stack Optimal Height under Uncertainty');

%% 
% There are many options available for the semi-infinite programming
% algorithm, |fseminf|. Consult the Optimization Toolbox(TM) User's Guide
% for details, in the Using Optimization Toolbox Solvers chapter, under
% Constrained Nonlinear Optimization: fseminf Problem Formulation and
% Algorithm.