Homework_05.ipynb

Kernel: Python 3 (system-wide)

Homework 05: Local Optimization

When asking questions about homework in Piazza please use a tag in the subject line like HW1.3 to refer to Homework 1, Question 3. So the subject line might be HW1.3 question. Note there are no spaces in "HW1.3". This really helps keep Piazza easily searchable for everyone!

For full credit, all code in this notebook must be both executed in this notebook and copied to the Canvas quiz where indicated.

Hints for 1-3

The code for this problem is very similar to the Wyndor Quadratic Program. The harder part is formulating the objective function and constraints.

Questions 1-3

The management of the Albert Hanson Company is trying to determine the best product mix for two new products. Because these products would share the same production facilities, the total number of units produced of the two products combined cannot exceed two per hour. Because of uncertainty about how well these products will sell, the profit from producing each product provides decreasing marginal returns as the production rate is increased. In particular, with a production rate of $R_1$ units per hour, it is estimated that Product 1 would provide a profit (in dollars per hour) of $200 R_1 - 100 R_1^2.$ If the production rate of product 2 is $R_2$ units per hour, its estimated profit (in dollars per hour) would be $300 R_2 - 100 R_2^2.$

Question 1 (8 points, manual)

Formulate a quadratic programming model and solve it using Pyomo with the ipopt solver. Put your complete code below and in your CoCalc notebook. We are allowing fractional solutions here, so this is not an integer programming problem.

Question 2 (1 point, auto)

Enter the maximum value of the profit to two decimal places.

Question 3 (1 point, auto)

To achieve the maximum value of the profit function, how many units of product 2 should be introduced each hour? Enter your answer to two decimal places.

Questions 4-7

Consider the nonconvex profit function $$p(x) = 100x^6 - 1359x^5 + 6836 x^4 - 15670 x^3 + 15870 x^2 - 5095 x$$ with $0 \leq x \leq 5.$

Question 4 (2 points, manual)

Graph the function for $0 \leq x \leq 5$ in your notebook and determine how many local minima there are. Do not count endpoints. Put your complete code to produce the graph below and in your CoCalc notebook.

Make sure your graph uses the correct range of values for x.

# we've defined the function for you to help avoid errors
profit = lambda x:100*x**6 - 1359*x**5 + 6836*x**4 - 15670*x**3 + 15870*x**2 - 5095*x

Question 5 (1 point, auto)

How many local minima are there? Do not count endpoints.

Question 6 (1 point, auto)

What is the value of $p$ at the second-largest local maximum? Enter your answer to two decimal places.

Question 7 (4 points, manual)

Write a multistart procedure that starts from uniform randomly sampled points in [0,5] to locate the absolute maximum value of profit. The algorithm should stop after 20 iterations in which no improvement in the max value has been obtained. Include your code and the output from your code in the space below. Your code should print out the max value and the x-value where it occurs.

This is pretty similar to the Rastrigin example in the section "Multistart for Finding Global Optima".
However you don't need the dim variable since there is only a single variable x in this problem.

Use x_initial = np.random.uniform(0,5) to get the initial search point in each iteration.

Questions 8-10

To find the line of least squares fit of the form $\hat{y} = b_0 + b_1 x$ to fit data of the form $(x_1,y_1), (x_2,y_2),\ldots,(x_n,y_n)$ we minimize a loss function. The loss function is the sum of the squares residuals and only depends on $b_0$ and $b_1$ for fixed $xy$-data:

$$SS(b_0,b_1) = \sum_{i=1}^{n} \left( y_i - \hat{y}_i \right)^2 = \sum_{i=1}^{n} \left( y_i - (b_0 + b_1 x_i) \right)^2$$

(Be sure to compute the squares before you compute the sum!)

The file age_height.csv contains ages (years) and heights (inches) for 7 children. Write Python code to evaluate the loss function (follow along with the logistic regression example while making suitable changes to the loss function) and use minimize to identify the coefficients of the line of least-squares fit for predicting height ($y$) from ($age$). Include a scatter plot of the data that includes a plot of the line.

Question 8 (10 points, manual)

Include complete Python code in the space below that includes the loss function, shows how to minimize it, prints out the output, and makes a plot of the line on a scatterplot of the data. You should also include your plot (check out savefig or use a cropped screenshot). The same elements should be in your CoCalc notebook. Here is the loss function from the logistic regression problem in the lesson along with a couple of notes to get you started.

def neg_log_loss( coef, *args): # rename your function
    b0 = coef[0]
    b1 = coef[1]
    x = args[0]
    y = args[1]
    
    # you'll need to change the code below here
    p = 1.0/(1.0 + np.exp(-(b0 + b1*x)))
    ll = sum( y*np.log(p)+(1-y)*np.log(1-p) )
    return(-ll) # log_loss needs to be MAXimized

Instead of p and ll compute ss the sum of the squares as shown in the formula above and then return(ss). Make sure you're computing the sum of the squared residuals and not the sum of the residuals squared ... parentheses matter!

If you need help plotting a line on a scatter plot you can modify the code in the lesson where you replace the sigmoid funciton with a function for your line.

This is also a perfect use case for ChatGPT.
I entered this prompt 'given x and y lists and the equation of a line show how to make a scatterplot with the x and y lists and plot the line on top using python' and got a useful response.

Question 9 (1 point, auto)

What is the value of the intercept, $b_0$? Enter your answer to three decimal places.

Question 10 (1 point, auto)

What is the value of the slope, $b_1$? Enter your answer to three decimal places.

Questions 11-12

The knapsack problem is a classical combinatorial optimization problem that will be good for practicing with the ideas of discrete local search and multistart. Given a set of items, each with a weight and a value, determine which items to include in a collection so that the total weight is less than or equal to a given limit and the total value is as large as possible. In the 0-1 version of the knapsack problem, the decision variables are binary (or boolean) and represent whether to include each item in the collection. We'll start with 20 items. You need to determine the collection of items that maximizes the value and keeps the total weight up to 50 (that is $\leq 50$).

# generate random weights and values for a knapsack problem (DO NOT CHANGE)
import numpy as np
num_items = 20
np.random.seed(seed=123)
values = np.random.randint(low=5, high=50, size=num_items)
weights = np.random.randint(low=1, high=10, size=num_items)
max_weight = 50
np.random.seed() # use system clock to reset the seed so future random numbers will appear random

The variables will be a vector of booleans of length num_items. We could initialize a vector like this and then set the vector to include the 1st, 3rd, and 5th items:

import numpy as np
x = np.zeros(num_items, dtype = bool)  # all false

array([ True, False,  True, False,  True, False, False, False, False,
       False, False, False, False, False, False, False, False, False,
       False, False])

The total weight of the items included in the collection:

tot_weight = sum( weights[x] )
tot_weight

11

The total value of the items included in the collection:

tot_value = sum( values[x] )
tot_value

68

Implement a local search where the search starts with no items included in the collection and generates new states (moves) by randomly choosing one of the booleans in the state vector and toggling it. Like this:

# try executing this cell a few times and watch x change
bit_to_flip = np.random.randint(num_items)
x[bit_to_flip] = ~x[bit_to_flip]
x

array([False, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False, False,  True, False,
       False, False])

Accept the state if the total weight is is $\leq 50$ and maximize the value by moving uphill

Question 11 (4 points, manual)

Write a local search algorithm that randomly changes one boolean and accepts the move if it improves the total value. Iterate until no improvements have been made in the last 1000 iterations. Write the algorithm as a function with the values and weights as inputs and that returns the best collection of items to include as well as the value and weight of that collection. Include your code and the results of running the code once in the space below.

The function you need to write has a similar structure to the random_reversal_search function in the lesson.
You'll have to adapt it work with the knapsack problem.

Question 12 (4 points, manual)

Implement a variation of multistart search from the lesson. Instead of starting from random initial states, start each local search from an empty knapsack. The random moves should include sufficient randomness. You should do 200 local searches. Each local search should call the function you defined in Question 11. Clearly identify the best overall solution. Put complete code to do the searches and display the best results in the space below and in your CoCalc notebook. Do not print out information from each search, just the final results.

Next week we'll see some alternative search techniques that will generally enable us to find better solutions.

Question 13 (2 points, manual)

Enter the maximum total value for knapsack that you found in Question 12. If it's not over 400 then something is likely wrong.

Question 14 (10 points, manual)

The Value Balancing Problem was a Self_Assessment problem near the end of the Lesson_05 notebook. You can see our solution in the Self_Assessement_Soln notebook. Use the locsearch package in this directory to find a local solution to the Value Balancing problem. You can follow the example at the end of Lesson_05. We've copied some code from the lesson to get you started. Remember to set this up so that your decision variables (the group assignments) are stored in the state attribute of your subclass.

# the objective function
def group_fitness(groups,num_groups,values):
    # groups must be a numpy array for this to work
    sums = np.array([ sum( values[ groups == i] ) for i in range(num_groups) ])
    return max(sums)-min(sums)

# the move function
def change_group(groups, num_groups, debug=False):
    #get a copy of the groups
    new_groups = groups.copy()    
    #select item to change
    switch = np.random.randint(0, groups.shape[0])
    #select new group value
    new_group = np.random.randint(0,num_groups)
    while groups[switch] == new_group:
        new_group = np.random.randint(0,num_groups)
    new_groups[switch] = new_group
    if debug:
        print(f'The item at {switch} should change to {new_group}')
        print(f'The initial groups are: {groups} and the changed groups are {new_groups}')
       
    return new_groups
    
# problem data
num_groups = 4
num_items = 1000
np.random.seed(5)
values = np.random.randint(2,20,size=num_items)
np.random.seed()

In the space below, include your complete code to solve the Value Balancing problem using the LocalSearcher class. You should also include output from a run of the code. The same elements should be in your ÇoCalc notebook.

We provide a partial solution below to help you get started. Try to understand this code as well as you can since we'll do a similar problem next week.

# If you're new to object oriented programming, we've provided a partial solution here for you to adapt

from locsearch import LocalSearcher
import numpy as np
import pandas as pd

# Define your class which inherits from LocalSearcher class
class ValueBalancingProblem(LocalSearcher):

    # pass in state and problem parameters (fix the XXXs)
    def __init__(self, state, values, XXX, debug):
        self.values = values
        self.XXX = XXX
        self.debug = debug
        super(ValueBalancingProblem, self).__init__(state)  # important!
        
    def move(self):
        self.state = # call the function that makes a move, pass in self.state, self.values, etc.
        
    def objective(self):
        return # call function that evaluates the objective function, pass in self.state, etc.
    
num_groups = 4
num_items = 1000
debug = False

# generate problem data (don't change)
np.random.seed(5)
values = np.random.randint(2,20,size=num_items)
np.random.seed()

# create initial state
initial_groups = # use numpy randint to generate a random initial groups assignment

# create an object of the ValueBalancingProblem class
bestValues = ValueBalancingProblem( # pass everthing needed by your __init__ method to construct the object )

# tweak some local parameters to more closely observe the search
bestValues.max_no_improve= 50
bestValues.update_iter = 10

# call the local search method in our object to actually do the search (replace XXX by ???)
best_groups, best_fitness = bestValues.XXX()

print(f"Max difference between groups is {best_fitness}")