%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % (CHEM-E7225) ADVANCED PROCESS CONTROL D % #02 - Nonlinear Optimization with Opti Stack % % CasADi is a tool to quickly implement and efficiently solve dynamic % optimisation (optimal control) problems. You can also think of it as % a general toolkit for gradient-based optimisation % % CasADI is open-source, LGPL-licensed % Download from https://web.casadi.org (current version: 3.6.4) % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% close all; clear; clc addpath('~/Documents/MATLAB/casadi') % --> addpath('~/change_me_with_path_to/casadi-v3.6.4') import casadi.* %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % UNCONSTRAINED OPTIMIZATION % % Define the optimization problem % % min (10*y - x^2)^2 % % and solve for (x*,y*) from initial solution (-1, 1) % opti = casadi.Opti(); % Creates an Opti Stack structure x = opti.variable(); % Creates a scalar decision variable (x) y = opti.variable(); % Creates a scalar decision variable (y) opti.minimize( (10*y - x^2)^2 ); % Define the objective function to be minimized opti.solver('ipopt'); % Chooses a solver (IPOPT, qpOASES, ...) opti.set_initial([x y], [-1 1]); % Sets an initial solution (default: 0) opti.solve(); % Executes the solver x_sol = opti.value(x); % Retrieves the optimal solution of (x) y_sol = opti.value(y); % Retrieves the optimal solution of (y) %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % CONSTRAINED OPTIMIZATION % % Define the optimization problem % % min (10*y - x^2)^2 % s.t. x^2 + y^2 = 1 % x + y >= 1 % % and solve for (x*,y*) from initial solution (-1, 1) % opti = casadi.Opti(); % Creates an Opti Stack structure x = opti.variable(); % Creates a scalar decision variable (x) y = opti.variable(); % Creates a scalar decision variable (y) opti.minimize( (10*y - x^2)^2 ); % Define the objective function to be minimized opti.subject_to( x^2 + y^2 == 1 ); % Defines an equality constraint opti.subject_to( x + y >= 1 ); % Defines an inequality constraint opti.solver('ipopt'); % Chooses a solver (IPOPT, qpOASES, ...) opti.set_initial([x y], [-1 1]); % Sets an initial solution (default: 0) opti.solve(); % Executes the solver x_sol = opti.value(x); % Retrieves the optimal solution of (x) y_sol = opti.value(y); % Retrieves the optimal solution of (y) %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % CONSTRAINED OPTIMIZATION (MATRIX FORM) % % Define the optimization problem % % min c^T x + 1/2 x^T B x % s.t. L <= x <= U % % with c = [1], B = [5 2], L = [-10], U = [10] % [2] [2 10] [-10] [10] % % and solve for x* from initial solution x* = 0. % % Optimization data c = [1; 2]; B = [5 2; 2 10]; L = [-10; -10]; U = [ 10; 10]; % Optimization definition opti = casadi.Opti(); % Creates an Opti Stack structure x = opti.variable(2,1); % Creates a (2 x 1) decision vector (x) opti.minimize( c'*x + 1/2 * x'*B*x ); % Defines the quadratic objective function opti.subject_to( L <= x <= U ); % Defines the bounds on the decision variables opti.solver('ipopt'); % Chooses a solver (IPOPT, qpOASES, ...) opti.set_initial(x, zeros(2,1)); % Sets an initial solution (default: zeros(n,1)) opti.solve(); % Executes the solver x_sol = opti.value(x); % Retrieves the optimal solution of (x) %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % CONSTRAINED OPTIMIZATION (SUMMATIONS) % % Define the optimization problem (note: a Linear Regression) % % min \sum_{i=1}^{N} (theta^T*x_i - y_i)^2 % s.t. \sum_{i=1}^{2} |theta_i| <= 1 % % with x_1,...,x_N ~ Normal(0,I_2) and y_1,...,y_N ~ Normal(0,1) % and solve for theta* (\in R^2). % % Optimization data N = 100; % Number of "datapoints" x = randn(2,N); % Datapoints: Inputs x_i \in R^2 y = randn(1,N); % Datapoints: Outputs y_i \in R % Optimization definition opti = casadi.Opti(); % Creates an Opti Stack structure theta = opti.variable(2); % Creates a (2 x 1) decision vector (x) J = 0; % Summation for cost (J) for i = 1:N % Iterates over all i \in [1,N] J = J + ( (theta'*x(:,i) - y(i))^2 ); % Adds the current loss |theta^T*x_i - y_i| to J(.) end % -- G = 0; % Summation for constraint (G) for i = 1:2 % Iterates over all i \in [1,2] G = G + ( abs(theta(i)) ); % Adds the current term abs(theta_i) to G(.) end % -- opti.minimize( J ); % Defines the quadratic objective function opti.subject_to( G <= 1 ); % Defines the bounds on the decision variables opti.solver('ipopt'); % Chooses a solver (IPOPT, qpOASES, ...) opti.solve(); % Executes the solver theta_sol = opti.value(theta); % Retrieves the optimal solution of (theta) %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % CONSTRAINED OPTIMIZATION (SUMMATIONS + TRACKING ITERATIONS) % % Alternative to above problem using vectorized operations % and tracking the intermediate solutions from the IPOPT iterations % % Optimization data N = 100; % Number of "datapoints" x = randn(2,N); % Datapoints: Inputs x_i \in R^2 y = randn(1,N); % Datapoints: Outputs y_i \in R % Optimization definition opti = casadi.Opti(); % Creates an Opti Stack structure theta = opti.variable(2); % Creates a (2 x 1) decision vector (x) J = sum( ( theta'*x - y ).^2 ); % Vectorized cost function G = sum( abs( theta ) ); % Vectorized constraint function opti.minimize( J ); % Defines the quadratic objective function opti.subject_to( G <= 1 ); % Defines the bounds on the decision variables % Adds a callback to 'track' the value of (x) theta_iter = []; opti.callback(@(i) evalin('base', 'theta_iter = [theta_iter opti.debug.value(theta)];') ) % - opti.solver('ipopt'); % Chooses a solver (IPOPT, qpOASES, ...) opti.set_initial(theta, [-0.5, 0]); % Sets an initial solution (default: zeros(n,1)) opti.solve(); % Executes the solver theta_sol = opti.value(theta); % Retrieves the optimal solution of (theta) % The function below will plot the cost function heatmap (masking out the % unfeasible set), together with the value of (theta) at each iteration of % IPOPT. You can use it as a template for Homework 2. plotOptimization(J, G-1, theta, theta_iter)