why don't we have a a good result from this model

Answered

Zhenting Xu

• September 27, 2024 20:12

u_ini = sdpvar(obj.M*obj.T_ini,1,'full'); % u_ini is the initial input

y_ini = sdpvar(obj.P*obj.T_ini,1,'full'); % y_ini is the initial output

u = sdpvar(obj.M*obj.horizon,1,'full');

y = sdpvar(obj.P*obj.horizon,1,'full');

g = sdpvar(obj.T-obj.T_ini-obj.horizon+1,1,'full');

slack_u = sdpvar(obj.T_ini*obj.P,1,'full');

slack_y = sdpvar(obj.T_ini*obj.M,1,'full');

U_p = obj.Up;

U_f = obj.Uf;

Y_p = obj.Yp;

Y_f = obj.Yf;

 

 

A = [U_p;Y_p;U_f;Y_f];

b = [u_ini+slack_u;y_ini+slack_y;u;y];

constraints = [];

constraints = [constraints,A*g == b];

obj.optimization_problem = OptimizationProblem(u_ini,y_ini,u,y,g,...

slack_u,slack_y,constraints,problem_loss);

 

u_ini = reshape(data_ini.u(1:obj.T_ini,:)',[],1);% initial Input

y_ini = reshape(data_ini.y(1:obj.T_ini,:)',[],1);% initial Output

assign(obj.optimization_problem.variables.u_ini, u_ini);

assign(obj.optimization_problem.variables.y_ini, y_ini);

 

options = sdpsettings('verbose', 1, 'solver', 'gurobi','debug',1,'gurobi.QCPDual', 1 ...

,'gurobi.TuneTimeLimit', 0);

g_value = value(obj.optimization_problem.variables.g);

Related to

• Performance

0

• 

  Ryuta Tamura

  Gurobi Staff

  • September 28, 2024 15:09
  • Edited

  Hi,

  Could you clarify what you mean by “good result”? There are several possible candidates: results of variables are not what you want, optimization is slow, etc. Also, can you provide details of the model you are trying to solve, execution logs, etc.? They help us understand the situation.

  Thanks
  Ryuta

   

  0
• 

  Zhenting Xu

  • September 30, 2024 20:42

  

  

  function obj = build_problem(obj,...

  lambda_g,lambda_y,lambda_u,lambda_proj)

  if nargin < 2, lambda_g = 0 ; end

  if nargin < 3, lambda_y = 0 ; end

  if nargin < 4, lambda_u = 0 ; end

  if nargin < 5, lambda_proj = 0 ;end

  % Assertions

  if exist('build_loss', 'file') == 2

  % nothing to do

  else

  disp('Loss function callback cannot be none');

  end

  assert(lambda_g >= 0 && lambda_y >= 0,'Regularizers must be non-negative');

  assert(lambda_u >= 0,'Regularizer of u_ini must be non-negative');

  assert(lambda_proj >= 0,'The projection regularizer must be non-negative');

  obj.optimization_problem = []; % optimization problem is not build yet

  %Build variables

  %'full' indicates that the type of the variable is a full asymmetric matrix

  u_ini = sdpvar(obj.M*obj.T_ini,1,'full'); % u_ini is the initial input

  y_ini = sdpvar(obj.P*obj.T_ini,1,'full'); % y_ini is the initial output

  u = sdpvar(obj.M*obj.horizon,1,'full');

  y = sdpvar(obj.P*obj.horizon,1,'full');

  g = sdpvar(obj.T-obj.T_ini-obj.horizon+1,1,'full');

  slack_u = sdpvar(obj.T_ini*obj.P,1,'full');

  slack_y = sdpvar(obj.T_ini*obj.M,1,'full');

  U_p = obj.Up;

  U_f = obj.Uf;

  Y_p = obj.Yp;

  Y_f = obj.Yf;

   

  if lambda_proj > DeePC.SMALL_NUMBER

  % Compute projection matrix (for the least square solution)

  Zp = [U_p;Y_p;U_f];

  ZpInv = pinv(Zp);

  I = eye(obj.T-obj.T_ini-obj.horizon+1);

  I_min_P = I-(ZpInv*Zp);

  end

   

  A = [U_p;Y_p;U_f;Y_f];

  b = [u_ini+slack_u;y_ini+slack_y;u;y];

  % Build constraints

  constraints = [];

  constraints = [constraints,A*g == b];

  if isclose(lambda_y,0)

  constraints = [constraints, norm(slack_y,2) <= DeePC.SMALL_NUMBER];

  end

  if isclose(lambda_u,0)

  constraints = [constraints, norm(slack_u,2) <= DeePC.SMALL_NUMBER];

  end

  %Reshape u and y

  %in order to calculate loss in all step k(1 to horizon)

  u = reshape(u,[obj.horizon, obj.M]);

  y = reshape(y,[obj.horizon, obj.P]);

  if exist('build_constraints','file')==2

  Constraints = build_constraints(u,y);

  else

  Constraints = [];

  end

  constraints = [constraints,Constraints];

  % Build Loss

  loss = build_loss(u,y);

  if lambda_g > DeePC.SMALL_NUMBER

  regulariers = lambda_g * norm(g,1); % add L1 norm of g

  else

  regulariers = 0;

  end

  if lambda_y > DeePC.SMALL_NUMBER

  regulariers = regulariers + lambda_y*norm(slack_y,1);

  end

  if lambda_proj > DeePC.SMALL_NUMBER

  % Add projection regularizer

  regulariers = regulariers + lambda_proj*norm(I_min_P*g);

  end

  if lambda_u > DeePC.SMALL_NUMBER

  % Add L1 norm of slack_u

  regulariers = regulariers + lambda_u*norm(slack_u,1);

  end

  % Define the total loss

  problem_loss = loss + regulariers;

  % Define the optimization problem

  obj.optimization_problem = OptimizationProblem(u_ini,y_ini,u,y,g,...

  slack_u,slack_y,constraints,problem_loss);

  end

  function [u_optimal,info] = solve(obj,data_ini,varargin)

  % Assertions

  assert(size(data_ini.u,2)==obj.M,'Incorrect number of features for the Input signal');

  assert(size(data_ini.y,2)==obj.P,'Incorrect number of features for the Output signal');

  assert(size(data_ini.u,1)==size(data_ini.y,1),'Input/Output data must have the same length');

  assert(size(data_ini.y,1)==obj.T_ini,'Invalid size');

  assert(~isempty(obj.optimization_problem),'Problem was not built');

  % Need to transpose to make sure that time is over the columns, and features over the rows

  u_ini = reshape(data_ini.u(1:obj.T_ini,:)',[],1);% initial Input

  y_ini = reshape(data_ini.y(1:obj.T_ini,:)',[],1);% initial Output

  assign(obj.optimization_problem.variables.u_ini, u_ini);

  assign(obj.optimization_problem.variables.y_ini, y_ini);

  % assign(obj.optimization_problem.variables.slack_u,0);

  % assign(obj.optimization_problem.variables.slack_y,0);

  options = sdpsettings('verbose', 1, 'solver', 'gurobi','debug',1,'gurobi.QCPDual', 1 ...

  ,'gurobi.TuneTimeLimit', 0);

  % options.mosek.MSK_DPAR_INTPNT_TOL_REL_GAP = 1e-8;

  % options.mosek.MSK_DPAR_INTPNT_TOL_PFEAS = 1e-8;

  % options.mosek.MSK_IPAR_INTPNT_MAX_ITERATIONS = 10000;

  try

  result = optimize(obj.optimization_problem.constraints, obj.optimization_problem.objective_function,options, varargin{:});

  catch ME

  error('Error while solving the DeePC problem.Details: %s',ME.message);

  end

   

  %If the problem is unbounded

  if result.problem == 1

  error('Problem is unbounded');

  end

   

  %value is a function provided by YALMIP

  %to extract the values of the optimized variables after the solution.

  g_value = value(obj.optimization_problem.variables.g);

  %compute the optimal input signal

  %In calculation:u1..u_horizon in one column, but reshape in

  %(u1...u_horiozon)^T

  u_optimal = reshape(obj.Uf*g_value,[obj.horizon,obj.M]);

         end

      end

  end

  0
• 

  Zhenting Xu

  • September 30, 2024 20:44

  

  

  function loss = build_loss(u,y)

    

      P = size(y, 2)

     y_ref = ones(size(y));

  % Q & R are the weight matrix

    loss = norm(y-y_ref, 2)^2;

  end

   

  function Constraints = build_constraints(u,y)

  horizon = size(u, 1);

  M = size(u, 2);

  P = size(y, 2);

  Constraints = [];

  % for i=1 : length(horizon)

  Constraints = [Constraints,-1<=u<=1];

  end 

  0
• 

  Zhenting Xu

  • September 30, 2024 20:44

  

  

  %% intialize the parameters in DeePC

  % DeePC Paramters

  s=1; % How many steps before we solve again the DeePC problem

  T_ini = 2;% Size of initial set of data

  %number of data points used to estimate the system

  T_list = [];

  T_list = [T_list,100];

  %Horizon length

  horizon = 10;

  LAMBDA_G_REGULARIZER = 0; % g regularizer

  LAMBDA_Y_REGULARIZER = 0; % y regularizer

  LAMBDA_U_REGULARIZER = 0; % u regularizer

  EXPERIMENT_HORIZON = 100; % Total number of steps

  A = [0.70469, 0;

  0.24664, 0.70469];

  B = [0.75937;

  0.12515];

  C = [0, 1];

  D = zeros(size(C, 1), size(B, 2));

  % state model aufstellen and sampling time is 1

  % sys is a class instance of systems

  sys = systems(ss(A, B, C, D, 1));

  %% Simulate for different values of T (Data set T times)

  for i = 1: length(T_list)

  fprintf('Simulating with %d initial samples...\n', T_list(i));

  sys=sys.reset();

  %Generate initial data and initialize DeePC

  [data,sys] = sys.apply_input(normrnd(0,1,[T_list(i),1]),0);

  deepc = DeePC(data, T_ini,horizon);

  %Create initial data (on the right side(constraint euqation))

  % u_ini & y_ini

  data_ini = Data(zeros(T_ini,1),zeros(T_ini,1));

  sys=sys.reset(data_ini);

   

  deepc = deepc.build_problem( ...

  LAMBDA_G_REGULARIZER, LAMBDA_Y_REGULARIZER, LAMBDA_U_REGULARIZER);

  for j = 1:(EXPERIMENT_HORIZON/s)

  % solve DeePC

  [u_optimal, info] = deepc.solve(data_ini);

  %Apply optimal control input

  [data,sys] = sys.apply_input(u_optimal(1:s,:),0);

  % Fetch last T_ini samples

  data_ini = sys.get_last_n_samples(T_ini);

  end

  data_all = sys.get_all_samples();

  figure;

  subplot(2,1,1);

  plot(data_all.y(T_ini:end), 'DisplayName', sprintf('s=%d, T=%d, T_i=%d, N=%d', s, T_list(i), T_ini, horizon));

  ylim([0, 1.5]);

  xlabel('t');

  ylabel('y');

  grid on;

  title('Closed loop - output signal y_t');

  legend('show', 'Location', 'best');

  subplot(2,1,2);

  plot(data_all.u(T_ini:end), 'DisplayName', sprintf('s=%d, T=%d, T_i=%d, N=%d', s, T_list(i), T_ini, horizon));

  ylim([-1.2, 1.2]);

  xlabel('t');

  ylabel('u');

  grid on;

  title('Closed loop - control signal u_t');

  legend('show', 'Location', 'best');

  sgtitle('Closed loop signals');

  end

  0
• 

  Zhenting Xu

  • September 30, 2024 20:46

  

  

  %% Define the new system

  classdef systems

  properties

  sys

  x0

  u

  y

  end

  methods

  function obj = systems(sys,x0)

  if nargin < 2

  x0 = zeros(1,size(sys.A, 1));

  end

  assert(~isempty(x0) && size(sys.A,1)==length(x0),'Invalid initial condition');

  obj.sys= sys;

  obj.x0 = x0;

  obj.u = [];

  obj.y = [];

  end

   

  function [data,obj] = apply_input(obj,u,noise_std)

  if nargin < 3

  noise_std = 0.5;

  end

  T = size(u,1);

  %How many times does it run u could be a matrix

  %The lsim function computes the response of the system given

  %input u and initial state obj.x0 and returns output y, time vector t, and state vector x0

  if T > 1

  [y_t,~,x_0] = lsim(obj.sys,u,(0:T-1)'*obj.sys.Ts,obj.x0);

  obj.x0 = x_0(end,:);% last column is final state as inital after T_ini

  else

  y_t = obj.sys.C*obj.x0';

  obj.x0 = (obj.sys.A*obj.x0' + obj.sys.B*u)';

  end

   

  y_t = y_t + noise_std*reshape(randn(T, 1), T, 1);

  if isempty(obj.u)

  obj.u = u;

  else

  obj.u = [obj.u; u];

  end

  % measurement noise on the output

  % measured value to the output

  if isempty(obj.y)

  obj.y = y_t;

  else

  obj.y = [obj.y; y_t];

  end

  data = Data(u,y_t);

  end

  function data = get_last_n_samples(obj,n)

  assert(size(obj.u,1) >= n,'Not enough samples are available');

  data = Data(obj.u(end-n+1:end,:),obj.y(end-n+1:end,:));

  end

  function data = get_all_samples(obj)

  data = Data(obj.u,obj.y);

  end

   

  %If data_ini or x0 are not provided, set them to empty arrays, respectively

  % here obj is also a paramter

  function obj = reset(obj,data_ini,x0)

  if nargin < 2

  data_ini = [];

  end

   

  if nargin < 3

  x0 = [];

  end

   

   

   

  if isempty(data_ini)

  obj.u = [];

  obj.y = [];

  else

  obj.u = data_ini.u;

  obj.y = data_ini.y;

  end

   

  if isempty(x0)

  obj.x0 = zeros(1,size(obj.sys.A,1));

  else

  obj.x0 = x0;

  end

  end

  end

  end

  Because the best control output u should be 1 at the beginning, but it always gives a result of 0.36

   

  0
• 

  Zhenting Xu

  • September 30, 2024 20:47

  

  

  %% intialize the parameters in DeePC

  % DeePC Paramters

  s=1; % How many steps before we solve again the DeePC problem

  T_ini = 2;% Size of initial set of data

  %number of data points used to estimate the system

  T_list = [];

  T_list = [T_list,100];

  %Horizon length

  horizon = 10;

  LAMBDA_G_REGULARIZER = 0; % g regularizer

  LAMBDA_Y_REGULARIZER = 0; % y regularizer

  LAMBDA_U_REGULARIZER = 0; % u regularizer

  EXPERIMENT_HORIZON = 100; % Total number of steps

  A = [0.70469, 0;

  0.24664, 0.70469];

  B = [0.75937;

  0.12515];

  C = [0, 1];

  D = zeros(size(C, 1), size(B, 2));

  % state model aufstellen and sampling time is 1

  % sys is a class instance of systems

  sys = systems(ss(A, B, C, D, 1));

  %% Simulate for different values of T (Data set T times)

  for i = 1: length(T_list)

  fprintf('Simulating with %d initial samples...\n', T_list(i));

  sys=sys.reset();

  %Generate initial data and initialize DeePC

  [data,sys] = sys.apply_input(normrnd(0,1,[T_list(i),1]),0);

  deepc = DeePC(data, T_ini,horizon);

  %Create initial data (on the right side(constraint euqation))

  % u_ini & y_ini

  data_ini = Data(zeros(T_ini,1),zeros(T_ini,1));

  sys=sys.reset(data_ini);

   

  deepc = deepc.build_problem( ...

  LAMBDA_G_REGULARIZER, LAMBDA_Y_REGULARIZER, LAMBDA_U_REGULARIZER);

  for j = 1:(EXPERIMENT_HORIZON/s)

  % solve DeePC

  [u_optimal, info] = deepc.solve(data_ini);

  %Apply optimal control input

  [data,sys] = sys.apply_input(u_optimal(1:s,:),0);

  % Fetch last T_ini samples

  data_ini = sys.get_last_n_samples(T_ini);

  end

  data_all = sys.get_all_samples();

  figure;

  subplot(2,1,1);

  plot(data_all.y(T_ini:end), 'DisplayName', sprintf('s=%d, T=%d, T_i=%d, N=%d', s, T_list(i), T_ini, horizon));

  ylim([0, 1.5]);

  xlabel('t');

  ylabel('y');

  grid on;

  title('Closed loop - output signal y_t');

  legend('show', 'Location', 'best');

  subplot(2,1,2);

  plot(data_all.u(T_ini:end), 'DisplayName', sprintf('s=%d, T=%d, T_i=%d, N=%d', s, T_list(i), T_ini, horizon));

  ylim([-1.2, 1.2]);

  xlabel('t');

  ylabel('u');

  grid on;

  title('Closed loop - control signal u_t');

  legend('show', 'Location', 'best');

  sgtitle('Closed loop signals');

  end

  0
• 

  Zhenting Xu

  • September 30, 2024 20:48

  

  

  function result = isclose(a, b,rel_tol,abs_tol)

  if nargin < 3

  rel_tol = 1e-9;

  end

  if nargin < 4

  abs_tol = 0.0;

  end

  result = abs(a - b) <= max(rel_tol * max(abs(a), abs(b)), abs_tol);

  end

   

  function [Up, Uf, Yp, Yf] = split_data(data, T_ini, horizon)

  % Utility function used to split the data into past data and future data.

  % Constructs the Hankel matrix for the input/output data, and uses the first

  % Tini rows to create the past data, and the last 'horizon' rows to create

  % the future data.

  %

  % :param data: A tuple of input/output data. Data should have shape TxM

  % where T is the batch size and M is the number of features

  % :param T_ini: number of samples needed to estimate initial conditions

  % :param horizon: horizon

  % :return: Returns Up,Uf,Yp,Yf

  assert(T_ini >= 1, 'Tini cannot be lower than 1');

  %check Tini whether >= 1

  assert(horizon >= 1, 'Horizon cannot be lower than 1');

  % check whether horizon >= 1

  [~,Mu] = size(data.u); %[u1;u2...] each row u has Mu dimensions

  % read the row & column numbers from the input u in the data

  [~,My] = size(data.y);

  % read .... from the output y

   

  %T_ini+horizon = oder data.u = col(u1...uT) data.y = col(y1...yT)

  Hu = create_Hankel_matrix(data.u,T_ini+horizon);

  Hy = create_Hankel_matrix(data.y,T_ini+horizon);

  % Past input daten Up

  Up = Hu(1:T_ini*Mu,:);

  % Input data in the future Uf

  Uf = Hu(end-horizon*Mu+1:end,:);

  % Past output data Yp

  Yp = Hy(1:T_ini*My,:);

  % output data in the future Yf

  Yf = Hy(end-horizon*My+1:end,:);

  end

  0
• 

  Zhenting Xu

  • September 30, 2024 20:50

  

  

  %% define the class

  classdef OptimizationProblemVariables

  properties

  u_ini

  y_ini

  u

  y

  g

  slack_u

  slack_y

  end

  methods

  function obj = OptimizationProblemVariables(u_ini,y_ini,u,y,g,slack_u,slack_y)

  obj.u_ini = u_ini;

  obj.y_ini = y_ini;

  obj.u = u;

  obj.g = g;

  obj.y = y;

  obj.slack_u = slack_u;

  obj.slack_y = slack_y;

  end

  end

  end

   

  classdef OptimizationProblem

  properties

  variables

  constraints

  objective_function

  end

  methods

  function obj = OptimizationProblem(u_ini,y_ini,u,y,g,slack_u,slack_y,constraints,objective_function)

  obj.variables = OptimizationProblemVariables(u_ini,y_ini,u,y,g,slack_u,slack_y);% call the constructor

  obj.constraints = constraints;

  obj.objective_function = objective_function;

  end

  end

  end

   

  classdef Data

  properties

  u;

  y;

  end

  methods

  function obj = Data(u,y)

  obj.u = u;

  obj.y = y;

  end

  end

  end

  0
• 

  Ryuta Tamura

  Gurobi Staff

  • October 01, 2024 00:45
  • Edited

  Hi,
  Could you give us the answer to the first question?

  Could you clarify what you mean by “good result”?

  What kind of good results are you looking for and what kind of results are you actually getting?
  This article may be helpful: How do I diagnose a wrong result? – Gurobi Help Center

  Thanks, 
  Ryuta

  0
• 

  Zhenting Xu

  • October 01, 2024 07:19

  

  

  My solution is wrong, the correct solution should be u=1

  thank you 

  0
• 

  Ryuta Tamura

  Gurobi Staff

  • October 01, 2024 08:33
  • Edited

  I found your another post seems to be related with this. You mentioned about "assign" there. According to the YALMIP reference, it seems to introduce initial values for variables.

  Note, assign does not replace or constrain the variable with a particular value. It is only used for initial values in nonlinear solvers.

  Therefore, it does not always lead to u=1. Are there any other factors leading to u=1?

  Thanks, 
  Ryuta

  0
• 

  Zhenting Xu

  • October 01, 2024 14:14

  

  

  So How to change the code 

  thank you very much            this initial value should be fixed until to next iteration

  0
• 

  Ryuta Tamura

  Gurobi Staff

  • October 01, 2024 16:09

  Hi,
  
  Does it not work if you write it as a normal constraint?
  
  Thanks, 
  Ryuta

  0
• 

  Zhenting Xu

  • October 01, 2024 20:30

  

  

  now  it's very useful  thank you very much    thanks 

   

  0

Please sign in to leave a comment.