Optimization Theory

If you have not previously used the Sympy package for symbolic programming, you should study the Sympy Intro module. 


The optimization theory that I implement in Python is for the constrained maximization problem with an equality constraint:  max F subject to G = c where F is the objective function to be maximized, G is the constraint function, and c is the constraint constant. 


I provide Python implementations of the tangency condition, Lagrange method, plotting the functions and level curves, and checking the functions for the quasi-properties of quasiconcavity and quasiconvexity.  To learn about optimization theory, I recommend Dixit (1970) and Morgan (2015).


Python Code


To receive the Python code for this topic, use the Contact form to request the zip file and I will e-mail it to you.


The optimization theory modules use the Sympy, Numpy, and Matplotlib packages as shown in the table below.






Dixit, A. (1990).  Optimization in Economic Theory. Oxford University Press.


Morgan, P. (2015).  An Explanation of Constrained Optimization for Economists.  University of Toronto Press.


Sympy Intro

Tangency Condition

Equality Constraints



Sympy Intro



In this Python module, you will learn some of the basic Python syntax and some Sympy methods and functions that are used in the optimization modules.  


New Python


Start:  The first task is to import packages that are used in the module and provide them with a reference name.  In this program we use only Sympy and refer to it using “sym”.  Any Sympy function is then called by prefacing the function name with “sym.” (sym period).  Sympy methods that operate directly on objects are not prefaced by “sym.” (sym period).


Case Sensitive:  It is important to realize that, like most computer languages, Python and all its packages are case sensitive.  For example, do not type X for x, or Symbols for symbols, and so on.


Section 1:  For symbolic computation, you need to tell Sympy what the symbols are.


Section 2:  Python actions are triggered by methods and functions.  A method operates on an object, and you call it by typing the object name, a period, the method name, and the parameters in ( ), like this:  object.method().  In this section, you use the Sympy subs( ) method.


Section 3:  Here you use some Sympy functions that take objects as parameters:  expand( ), factor( ), and simplify( ).


Section 4:  You will do a lot of differentiation, so this section shows you how to use the diff( ) function.


Section 5:  You will use lists, tuples, and finitesets, which are three of the data structures for sequences. 


List:  Elements of a list are contained in square brackets [ ] and indexed 0, 1, 2, etc.  You get an item from a list using its sequential index, e.g., myList[0].  You can add, delete, or change items in a list after it is defined, i.e., it is mutable.  Note:  The first element in a list has an index of 0, not 1!


Tuple: Elements are contained in round brackets ( ) and indexed like a list.  However, a tuple is immutable, i.e., not changeable after it is defined.


FiniteSet:  Elements are in round brackets ( ).  FiniteSets, unlike lists or tuples, cannot have multiple occurrences of the same element.  Elements of a set can be any immutable type, so elements can be numbers, strings, or tuples but not lists. 


We use two Sympy functions for solving equations:


solveset( ):  to solve a single-variable equation

nonlinsolve( ):  to solve a system of n nonlinear equations with n variables


These two functions produce output in a FiniteSet data structure.  The contents of the FiniteSet, however, depend on the function.


solveset( ):  The output for solveset( ) is a FiniteSet that contains the solution values.  The easiest way to get the results from the FiniteSet is to convert to a list as shown in line 73 and then get the individual solutions by index from the list. 


nonlinsolve( ):  The output for nonlinsolve( ) is also a FiniteSet, but it contains tuples that in turn contain the solution values.   The easiest way to get the results from the FiniteSet is to use the for iterator as shown in the code in line 85.

Tangency Condition



The purpose of this module is to demonstrate Python code for the tangency condition.  


This simple optimization method can be used to solve a two-variable (X1, X2) optimization problem with an objective function F subject to a constraint G = c.  The tangency condition is F1/F2 = G1/G2, where F1 refers to the partial derivative of F with respect to X1, and so on.  Solving the tangency condition results in X2 as a function of X1, which defines a set of possible optimal points.  You find the optimal point in this set by using the constraint G = c.


New Python


This program uses much of the Python that you learned in the Sympy Intro module.  This module is structured into six steps.  There are two new functions and some new parameters introduced in steps 2, 3, and 6:


Step 2:  The “sep” parameter is used now in the print function to separate individual elements by a comma.


Step 3:  When you need an equality condition, you use the Sympy Eq() function.  For example, for x = y use Eq(x, y).


Step 6:  Create a tuple to represent the optimal solution xbar by simply putting the elements for x1bar and x2bar into round brackets.  Alternatively, you could use a list by putting the elements into square brackets.  Use the evalf() function to change fractions to floating point and to specify the number of significant digits.

Equality Constraints



This module demonstrates Python code for solving an optimization problem with equality constraints using the Lagrangean.  


New Python


This module is structured into five steps.  There is some new Python in steps 4  and 5:


Step 4:  You "could" define first order conditions as follows:

FOC1 = sym.Eq(L1, 0)

FOC2 = sym.Eq(L2, 0)

FOClam1 = sym.Eq(Llam1, 0)


However, recall from the Sympy Intro module that the Sympy solve function assumes you want to find the zeros of the function(s).  Thus, instead of using sym.nonlinsolve([FOC1, FOC2, FOClam1], [x1, x2, lam1]), you can simply use sym. nonlinsolve ([L1, L2, Llam1], [x1, x2, lam1]) which Sympy interprets as sym. nonlinsolve ([L1=0, L2=0, Llam1=0], [x1, x2, lam1]).


Step 5:  There are many solutions to this problem that are output in a FiniteSet of tuples.  We must iterate through the FiniteSet to find the nonnegative solution.   


We use the “for” iterator to go through the set, checking for the solution with values greater than or equal to zero.  We do this elegantly by simply checking for a solution with a minimum value greater than or equal to zero.


Once we find a non-negative solution, we “break” out of the iteration.




In this module, 3-dimensional and 2-dimensional plots are added to the Equality Constraint module. 


You need two packages to do plotting, Numpy and Matplotlib.Pyplot.  Add the two import statements after the Sympy import statement. 


New Python


This module is structured into nine steps (0 to 8).  There are new Python concepts and functions in steps 0, 1, 6, 7, and 8:


Step 0:  The code for finding the non-negative tuple in a FiniteSet is very different from what we are doing in the rest of the program.  For this type of code, it is best to define a function.


All functions should be defined at the beginning of the module, before the rest of the code.  A function begins with “def” (which means “define”), and ends with a return of the results of the function.   


Step 1:  For plotting the F and G functions we need to calculate the functions over a grid of points.  For this calculation, we need to create Numpy functions from the symbolic Sympy functions.  This type of single expression function is called a “lambda” function, and it is elegantly created by using the Sympy lambdify() function.  


Step 6:  In Numpy you will use the linspace and meshgrid functions to create the canvas upon which you will paint your plot.  Note that the program prints the meshgrid and the values of F on the meshgrid, providing insight into how the plots are created.  You use several Pyplot functions, which are explained in the code.


Steps 7 and 8:  When adding subplots in Step 7 and Step 8, the 111 refers to the first subplot in a 1 by 1 figure.  The first two digits represent the number of rows and columns and the third digit represents the index.  I prefer in this case to have the plots in separate figures, but you could could put both plots in a 1 by 2 figure by adding subplots 121 and 122.




Figure produced by Plotting module:


Blue:  Level curves for the objective function F. 

Red:   Level curve of the constraint G = c.


Green star:  The optimal solution.












         Quasiconcave function                  Convex level curves                Quasiconvex function                Concave level curves




This module will add some matrix functions to your knowledge of Python. 


The quasi-properties of quasiconcavity and quasiconvexity are used in economic optimization to assess whether a solution is global and/or unique.   Briefly, the quasi-properties are determined by the level curves.  What is useful to us is the geometry of a function’s level curves, not the function itself.  Are the level curves strictly concave or strictly convex?  If they are strictly concave, the function is strictly quasiconvex and if they are strictly convex the function is quasiconcave.


New Python


There are new Python concepts and functions in steps 1, 2, and 3 of this module:


Step 1:  Because testing quasi-properties involves testing the sign of determinants, we tell Sympy that the variables are positive.


Step 2:  Sympy provides many Matrix functions, and here we use functions to insert the border into the Hessian.  Note that I am leaving the gradient as a row vector.  First, we create a bordered gradient by inserting a “column” which is just zero.  Second, we add the gradient as the last row.  Third, we add the transpose of the bordered gradient as the last column.  Note that we could have reversed this order and added the transpose of the gradient as the last column and then the bordered gradient as the last row.


Step 3:  We use the minor_submatrix() function to remove the 0 row and 0 column to create the C1 submatrix, calculate its determinant using the det() function, and determine the sign using the is_positive function.

Steps 4-6:  I copied the plotting code from the Plotting module and edited it for this function.