当前位置：网站首页>[mathematical modeling - Summary of planning model] | matlab solution

This paper summarizes the commonly used mathematical programming models in mathematical modeling , With detailed MATLAB Solve cases . It is divided into four modules ：

The general steps to solve the mathematical model are as follows ：

• Reading questions + Understand the model ;

• Design The decision variables ;

• List Objective function And constraint condition ;

• Express （ draw ） The feasible region represented by the constraints ;

• Find the optimal solution and value of the objective function in the feasible region ;

If you don't understand , Follow the example ...

1 Linear programming problem （LP）

seek `Linear objective function` stay `Linear constraints` The problem of maximum or minimum value under , They are collectively called linear programming problems .

Take a simple case ：

example ： A machine tool factory produces a 、 B two kinds of machine tools , The profit after sales of each unit is 4000 Yuan and 3000 element . The production of machine tool a needs A、B Machining , The processing time is... For each machine 2 Hours and 1 Hours ; The production of machine tool B needs A、B、C Three kinds of machining , The processing time is one hour for each set . If the machine hours available for processing each day are A machine 10 Hours 、B machine 8 Hours and C machine 7 Hours , Q: the factory should produce a 、 B. how many machine tools are there , To maximize the total profit ？

Establish a mathematical model according to the above description ：

Set up this factory to produce x1 Taijia machine tool and x2 The total profit of machine tool B is the largest , be x1,x2 Should satisfy :

In the above formula ：

• x1 , x2 Decision variables are called decision variables ;
•（1） The equation is called the objective function of the problem ;
•（2） Several inequalities in are the constraints of the problem , Write it down as s.t.( namely subject to);

The above objective functions and constraints are linear functions , So it is called linear programming problem .

When solving practical problems , Boil the problem down to A mathematical model of linear programming It's very important , But it is often the most difficult , Whether the model is properly established , It directly affects the solution . Choose the right The decision variables , It is one of the keys to establish an effective model .

The objective function of linear programming can be to find the maximum , You can also find the minimum value , The inequality sign of constraints can be less than or greater than . In order to avoid the inconvenience caused by this diversity of forms ,MATLAB The standard form of linear programming is :

among c and x by n Virial vector , A 、 Aeq Is a matrix of proper dimension ,b 、beq Is a column vector of proper dimension ;lb And ub Respectively x Lower and upper bounds of , Dimension and x Corresponding .

The model has been expressed as MATLAB The standard form of , You can use MATLAB Solved , Here the linprog function .

style 1

function lp_% MATLAB Solving linear programming problems ,linprog% c = -[4000;3000]; %  Objective function coefficient A = [2 1;1 1;0 1]; b = [10;8;7];Aeq = [];beq = [];lb = [0 0]; %  Lower bound of variable x = linprog(c,A,b,Aeq,beq,lb)fval = -c'*x

It should be noted that the constraints are transformed into a matrix with appropriate dimensions ！

If you don't want to matrix （ coefficient ） transformation , You can also use another programming style ：

style 2

%  Back to fval It's negative ,%  Even if the problem is positive ,%  stay matlab Inside ,prob2struct The maximization problem is transformed into the minimization problem with negative objective function x1 = optimvar('x1','LowerBound',0,'UpperBound',inf);x2 = optimvar('x2','LowerBound',0,'UpperBound',7);prob = optimproblem('Objective',4000*x1 + 3000*x2,'ObjectiveSense','max');prob.Constraints.c1 = 2*x1 + x2 <= 10;prob.Constraints.c2 = x1 + x2 <= 8;problem = prob2struct(prob);[sol,fval,~,output] = linprog(problem)

I believe that through the above simple examples , You have a certain understanding of modeling and solving linear problems , Now let's look at a more complex model .

2 Nonlinear programming

If the objective function · or · The constraints contain nonlinear functions , This kind of programming problem is called nonlinear programming problem .

There are two kinds of nonlinear programming problems: unconstrained and constrained , Their solution ideas are not quite consistent , Specific methods can be referred to “ Operations research ” perhaps “ optimization ” Books ; One of the special problems is Second planning problem .

As a general rule , Solving nonlinear programming is much more difficult than solving linear programming . and , Unlike linear programming Simplex method This general method , At present, there is no general algorithm for nonlinear programming , Each method has a specific scope of application .

The difference between linear programming and nonlinear programming is ： If the optimal solution of linear programming exists , Then the optimal solution can only be reached on the boundary of its feasible region （ In particular, the vertex of the feasible region reaches ）, And the optimal solution of nonlinear programming （ If there is ） It is possible to reach . The graphical method of the lower linear problem can easily reflect the property that the solution is at the vertex ：

Graphic method

The nonlinear programming problem is MATLAB The standard form in is ：

The commonly used function in nonlinear solution is fmincon, Take an example model to see its usage :

Below is fmincon Call format for ：

 Does not contain nonlinear constraints [X,FVAL,EXITFLAG,OUTPUT,LAMBDA] = fmincon(FUN,X0,A,B,Aeq,Beq,LB,UB)

 Contains nonlinear constraints [X,FVAL,EXITFLAG,OUTPUT,LAMBDA] = fmincon(FUN,X0,A,B,Aeq,Beq,LB,UB,NONLCON)

 Define the objective function and nonlinear constraint by yourself [X,FVAL,EXITFLAG,OUTPUT,LAMBDA] = fmincon(@(X)MYOBJ(X),X0,A,B,Aeq,Beq,LB,UB,@(X)MYCON(X)) Sex constraint function %  Objective function function F = MYOBJ(X)F = ......%  Nonlinear constraint function function [G,Heq] = MYCON(X)G = ..... %  Nonlinear inequality constraints ,G<=0Heq = ..... %  Nonlinear equality constraints ,Heq == 0

（ The code of the model is downloaded at the end of the article ）

use fmicon The solution result of is ：

In addition to the default Interior point method Outside , It can also be used for fmincon Function sets different optimization algorithms ：

Next is the dynamic programming problem .

3 Dynamic programming

Basic idea of Dynamic Planning

• Divide the multi-stage decision-making process into stages , Choose state variables appropriately 、 Decision variables to define the optimal index function , Thus, the problem is transformed into a family of sub problems of the same type , Then solve them one by one ;

• The solution starts from the boundary conditions , Proceed in reverse order , Step by step recursive optimization ． When solving each sub problem , We should use the optimal result of the subproblem it has worked out before ． The optimal solution of the last subproblem , Is the optimal solution of the whole problem ;

• The dynamic programming method is to separate the current segment from the future segments , An optimization method that combines current benefits with future benefits , Therefore, the optimal decision-making selection of each segment is considered from the overall situation , It is generally different from the optimal choice of this paragraph ;

Learn dynamic programming , First of all, we need to understand the multi-stage decision-making problem . such as ：

Production decision problem ： Enterprises in the production process , Because demand changes over time , Therefore, in order to obtain the best production benefit of the whole year , It is necessary to determine the production plan month by month or quarter by quarter according to inventory and demand in the whole production process .

Machine load distribution problem ： A certain machine can produce under high and low loads . It is required to formulate a five-year plan , At the beginning of each year , Decide how to redistribute the quantity produced by the intact machine under two different loads , Make the total output of products reach the maximum within five years .

Space shuttle flight control problem ： Because the environment of space shuttle movement is constantly changing , Therefore, according to the situation of the space shuttle flying in different environments , Constantly determine the direction and speed of the space shuttle （ state ）, So that it can save the most fuel and complete the flight mission （ Such as soft landing ）.

　　 For the optimization problem of multi-stage decision-making process , American mathematician Bellman Et al. 20 century 50 In the early s, the famous optimization principle was put forward , The multi-stage decision-making problem is transformed into a series of single-stage optimization problems , So as to solve one by one , A new method to solve this kind of process optimization problem has been established ： Dynamic programming . For the best path （ The best decision-making process ） All the stages , The path from the beginning of each stage to the end of the whole process , It must be the best of all possible paths from the beginning of this stage to the end of the whole process （ Optimal decision ）, This is it. Bellman Put forward the famous The principle of optimization . namely ： The sub strategy of an optimal strategy must also be optimal .

Let's look at a solution case directly :

example ： Equipped with m=4 Goals , Target value （ Importance and harmfulness ） Each are not identical , Numerical value Ak（k=1,2,…,m) Express , Plan to use n=5 Missile raid , Probability of missile hitting the target
in a_k Is constant , It depends on the characteristics of the missile and the nature of the target ;u_k For the number of ground missiles fired at the target , How to make a plan to maximize the expected shock effect ？

First, we need to build a model of the problem ：

The function file is established according to the model and dynamic programming algorithm ：

 The main function clear n = 5;x1 = [n;nan*ones(n,1)];x2 = [0:n]';x = [x1 x2 x2 x2][p f] = dynprog(x,'DecisF1','SubObjF1','TransF1','ObjF1')disp('p The third column of is the value of the decision variable ')

 Decision function function u = DecisF1(k,x)%  Decision function if k ==4    u =x;else    u = 0:x;endend

 State transition equation function y = TransF1(k,x,u)%  State transition equation y = x-u;end

 Stage k Index function of function v = SubObjF1(k,x,u)%  Stage k Index function of a = [0.2 0.3 0.5 0.9];A = [8 7 6 3];v = A(k)*(1-exp(-a(k)*u));v = -v;end

 Functions in basic equations function y = ObjF1(v,f)%  Functions in basic equations y = v+f;% y = -y;end

 Dynamic programming dynprog function function [p_opt, fval] = dynprog(x,DecisFun, SubObjFun, TransFun, ObjFun)%  Dynamic programming reverse order algorithm % x It's a state variable , A column represents a stage state ;% DecisFun(k,x) By stage k State variable of x Find out the corresponding allowable decision variables ;% SubObjFun(k,x,u) Is the stage index function ;% TransFun(k,x,u) Is the state transition function , among x It's the stage k Some state variable of ,u Is the corresponding decision variable ;% ObjFun(v,f) It's No k Index function from stage to final stage , When ObjFun(v,f) = v+f when , Input ObjFun It can be omitted ;% fval It's a column vector , Each element represents p_opt Corresponding initial state of each optimal policy group x The optimal function value of

In dynamic programming problems , Path planning And very important , Below is a list about A Star algorithm and Probability roadmap （ be based on dijkstra Algorithm ） Solution case of :

A Star algorithm

function path=AStar(obstacle,map)path=[];        % For storing paths open=[];        %OpenList, Nodes considered close=[];       %CloseList, Nodes not considered findFlag=false; %findFlag Used to judge while Whether the cycle is over %% 1. Put the starting point at Openlist in % open：[ Node coordinates , On behalf of value F=G+H, On behalf of value G, Parent node coordinates ]open =[map.start(1), map.start(2) ,...    0+h(map.start,map.goal) ,...    0 , ...    map.start(1) , map.start(2)];% The coordinate difference with the current node （ The first two columns ）% The distance from the current node （ The last column ）NEXT = [-1,1,14;...    0,1,10;...    1,1,14;...    -1,0,10;...    1,0,10;...    -1,-1,14;...    0,-1,10;...    1,-1,14];%% 2.  Cycle through the following process while ~findFlag         %-------------------- First, judge whether the target point is reached , Or no path -----    if isempty(open(:,1))        disp(' There is no target path ');        return;    end        % Determine whether the target point appears in open In the list     [isopenFlag,Id]=isopen(map.goal,open);    if isopenFlag        disp(' Find the target ');        close = [open(Id,:);close];        findFlag=true;        break;    end            % according to Openlist The third column in （ Cost function F） Sort , lookup F The node with the lowest value     [Y,I] = sort(open(:,3)); % Yes OpenList Sort the third column in     open=open(I,:);%open The node in the first row in the F Minimum value         % take F The node with the lowest value ( namely open In the first line of the node ), Put it in close first line (close It's a constant backlog ), As the current node     close = [open(1,:);close];    current = open(1,:);    open(1,:)=[];% Because it's already from open Removed from , So the first column needs to be empty         % For the current node around 8 Calculate adjacent nodes , The main body of the algorithm ：------------------------    for in=1:length(NEXT(:,1))                        m=[current(1,1)+NEXT(in,1) , current(1,2)+NEXT(in,2) , 0 , 0 , 0 ,0];        m(4)=current(1,4)+NEXT(in,3); % m(4)   Adjacent nodes G value         m(3)=m(4)+h(m(1:2),map.goal);% m(3)   Adjacent nodes F value                         if isObstacle(m,obstacle)            continue;        end        [flag,targetInd]=findIndex(m,open,close);        if flag==1            continue;        elseif flag==2            m(5:6)=[current(1,1),current(1,2)];            open = [open;m];        else            if m(3) < open(targetInd,3)                m(5:6)=[current(1,1),current(1,2)];                open(targetInd,:) = m;            end        end    end        % draw     pause(0.01);    fillPlot(close,[204,159,42]/255);    hold on;    fillPlot(open,[51,102,153]/255);    end%% 3.  Tracing path path=GetPath(close,map.start);end

be based on dijkstra Probability roadmap for

%  Load map % --------------------------------------------------------—————————% PIC = imread(uigetfile('*.png'));         %  Load map PIC = imread('map2.png');                    %  Load map imshow(PIC);                                %  Show map hold on;%  The key parameters % --------------------------------------------------------—————————feasibleR = 51;                             %  The red value of the feasible area 【0~255】NUM = 200;                                  %  Discrete point scale %  Some style settings , You can comment out % --------------------------------------------------------—————————ax = gca;ax.Parent.Units = 'normalized';ax.Parent.Color = 'k';ax.Parent.Position = [0.1 0.1 0.8 0.7];ax.Units = 'normalized';ax.Position = [0.1 0.001 0.8 0.998];ax.Color = [51 51 51]/255;%  Select from / Starting point % --------------------------------------------------------—————————StartPoint = ginput(1);EndPoint = ginput(1);                       scatter([StartPoint(1) EndPoint(1)],[StartPoint(2) EndPoint(2)],...    'filled','SizeData',150,'MarkerFaceColor',[199 237 53]/255,'MarkerEdgeColor','w');%  Draw the starting point pause(1)%  Generate and draw feasible points %__________________________________________________________________________Height = size(PIC,1);Width = size(PIC,2);randPoints = floor([rand(NUM,1)*Height,...    rand(NUM,1)*Width])+1;                  %  Scatter points randomly in space allPoints = [StartPoint;EndPoint];          %  Initialize the feasible point list for i=1:NUM    if(PIC(randPoints(i,2),randPoints(i,1),1) == feasibleR  )        allPoints = [allPoints;double(randPoints(i,:))];    endend%  Get a list of all possible points for i = 1:size(allPoints,1)    scatter(allPoints(i,1),allPoints(i,2),'filled','SizeData',20,'MarkerFaceColor',[199 237 53]/255);    pause(.0005)end%  Initialize the adjacency matrix of undirected graph %__________________________________________________________________________AdjacencyMatrix = genAM(PIC,allPoints);%  The matrix AdjacencyMatrix Use 【Dijkstra Shortest path algorithm 】;%__________________________________________________________________________

4 Multiobjective programming

The above problem, no matter how complex the constraints are , There is only one goal . But in reality, many problems can finally be abstracted as Multi objective optimization problem . This is because practical problems are often complicated , Generally, many aspects will be weighed when reaching the goal . Here is a brief introduction to several mainstream methods .

Multi objective problem model

• Contains multiple objective functions that may conflict • Hope to find a solution set that can well balance all optimization objectives

Understand multi-objective optimization , You need to understand first Pareto is the best （Pareto Optimal）

Pareto is the best

Pareto optimality refers to an ideal state of resource allocation . Given an inherent group of people and allocable resources , If there is a change from one allocation state to another , Without making anyone worse , Make at least one person better , This is Pareto improvement .

Pareto's optimal state is that there can be no more Pareto improved state ; let me put it another way , It is impossible to improve the situation of some people without hurting anyone else .

control （Dominace）

When x1 and x2 When the following conditions are met, it is called x1 control x2： For all objective functions x1 No comparison x2 Bad ; On at least one objective function ,x1 Strict ratio x2 It is better to ;

( The midpoint in the figure above 1 Dominant point 2; spot 5 Dominant point 1; spot 1 Sum point 4 They don't dominate each other )

Non dominated solution set

（Non-dominated solution set） When any solution in a solution set cannot be dominated by other solutions in the set , Then the solution set is called non dominated solution set .

Pareto optimal solution set

（Pareto-optimal set） The non dominated solution set of all feasible solutions is called Pareto optimal solution set .

Pareto optimal frontier

（Pareto-optimal front） The boundary of Pareto optimal solution set （boundary） It is called Pareto optimal frontier .

What are the classical methods for multi-objective optimization problems ？

Linear weighting method

The weight represents the importance of each objective function . The advantage is simplicity ; The disadvantage is that it is difficult to set a weight vector to obtain the Pareto optimal solution ; In some non convex cases, the Pareto optimal solution cannot be guaranteed .

Constraint transformation method

Multi objective genetic algorithm

Multi-objective optimization based on genetic algorithm is to use the principle of genetic algorithm to search the Pareto optimal frontier . The introduction of genetic algorithm has been pushed in previous periods , Search history push directly .

The basic process is as follows ：

• Randomly generate the initial population ;

• Calculate the objective function value and constraint function value of each point ;

• Grade the group according to the value of the objective function ;

• Calculate the constraint penalty of each point according to the value of the constraint function and the grading result 、 Bad penalty and total penalty ;

• Calculate the fitness according to the total penalty of each point ;• According to the fitness of each point , Make a selection 、 Crossover and mutation operations , Generate new groups ;

• The total penalty is 0 The points of are put into the non inferior solution set candidate table , Check the candidate list , Keep the second 1 Grade is not inferior , Delete other points ;

• Check for convergence , If not , Go to step 2;

• Delete the points in the candidate table that are too close to other points ;

• Output Pareto Optimal solution set and corresponding objective function value ;

• The decision-maker starts from Pareto The best solution set selects the most suitable solution for the problem ;

The advantage of genetic algorithm compared with traditional algorithm is that it can get an optimal solution set , Not just an optimal solution , This gives us more choices . But the computational complexity may be slightly higher , And some functions involved in it need careful design .

原网站

版权声明
本文为[Love coldplay]所创，转载请带上原文链接，感谢
https://yzsam.com/2022/207/202207260314007428.html