refactor: Restructure Repository to add eta optimization

Signed-off-by: Tuan-Dat Tran <tuan-dat.tran@tudattr.dev>

00_aoi_caching_simulation/note.md

@@ -83,4 +83,15 @@ CPU times: user 3min 46s, sys: 43 s, total: 4min 29s
 Wall time: 4min 29s
 for ACCESS_COUNT_LIMIT = 10_000 # Total time to run the simulation
 
-## Notes 11/27/2024
+## Notes 11/29/2024
+
+C_m = cost for cache miss
+C_delta = cost for refresh
+C = Cm + C_delta over all objects summarized
+
+We wanna minimize cost function
+
+N = number of objects
+B is cache size
+
+C_f roughly equals C_m

01_nb_cncf_optimization/.ipynb_checkpoints/nb_cost_optimization-checkpoint.ipynb

@@ -0,0 +1,6 @@
+{
+ "cells": [],
+ "metadata": {},
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

01_nb_cncf_optimization/matlab/.ipynb_checkpoints/Main-checkpoint.ipynb

@@ -0,0 +1,258 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "id": "d0996120-bb17-4476-b912-ce155100b2cb",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Hit probabilities (numerical and theoretical):\n",
+      " [[0.4  0.  ]\n",
+      " [0.4  0.08]\n",
+      " [0.4  0.16]\n",
+      " [0.4  0.24]\n",
+      " [0.4  0.32]\n",
+      " [0.4  0.4 ]\n",
+      " [0.4  0.48]\n",
+      " [0.4  0.56]\n",
+      " [0.4  0.64]\n",
+      " [0.4  0.72]\n",
+      " [0.4  0.8 ]]\n",
+      "Objective function values (numerical, theoretical): [33.17, 22.944080000000007]\n",
+      "Constraint violations (numerical, theoretical): [0.0, 0.0]\n"
+     ]
+    }
+   ],
+   "source": [
+    "import numpy as np\n",
+    "from scipy.optimize import minimize\n",
+    "import numpy as np\n",
+    "\n",
+    "# Define Parameters\n",
+    "lambda_vals = np.array([0.03, 0.04, 0.05, 0.06, 0.07, 1, 1.1, 1.2, 1.3, 1.4, 1.5])  # Request rates ascendingly\n",
+    "N = len(lambda_vals)\n",
+    "B = 4.4  # Cache size\n",
+    "c_delta = 1  # Age linear cost\n",
+    "c_f = 7  # Fetching linear cost (cache miss cost)\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "fcf0c13c-5b2c-457e-9aa6-8d349fcf13fa",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "\n",
+    "def theoretical_opt(lambda_vals, B, c_f, c_delta):\n",
+    "    \"\"\"\n",
+    "    Perform theoretical optimization to compute optimal hit probabilities.\n",
+    "    \n",
+    "    Parameters:\n",
+    "    - lambda_vals: Array of request rates.\n",
+    "    - B: Cache size (constraint on total hit probabilities).\n",
+    "    - c_f: Cost of fetching (cache miss cost).\n",
+    "    - c_delta: Cost of caching (hit cost).\n",
+    "\n",
+    "    Returns:\n",
+    "    - h_optt: Optimal hit probabilities.\n",
+    "    \"\"\"\n",
+    "    N = len(lambda_vals)\n",
+    "    h_optt = np.zeros(N)  # Initialize optimal hit probabilities\n",
+    "    differenc_set = np.arange(N)  # Set of variables to optimize\n",
+    "    fix_i = []  # Set of fixed variables (those already optimized)\n",
+    "    n = N\n",
+    "    b = B\n",
+    "    flag = True\n",
+    "\n",
+    "    while flag:\n",
+    "        if n == 0:  # If no variables left to optimize\n",
+    "            if b > 0:  # If there is leftover cache size, redistribute it\n",
+    "                differenc_set = np.where(h_optt == 0)[0]  # Find zero hit probability variables\n",
+    "                fix_i = np.setdiff1d(np.arange(N), differenc_set)\n",
+    "                n = len(differenc_set)\n",
+    "                continue\n",
+    "            else:  # No variables to optimize, finalize\n",
+    "                h_optt[differenc_set] = 0\n",
+    "                break\n",
+    "        \n",
+    "        # Calculate the optimal Lagrange multiplier (mu) and hit probabilities for the set of variables\n",
+    "        mu = max(0, (n * c_f - b * c_delta) / np.sum(1 / lambda_vals[differenc_set]))\n",
+    "        h_optt[differenc_set] = (c_f - mu / lambda_vals[differenc_set]) / c_delta\n",
+    "        \n",
+    "        # If mu < 0, adjust the cache size to set mu to zero in the next iteration\n",
+    "        if mu < 0:\n",
+    "            b = (n * c_f / c_delta)\n",
+    "            continue\n",
+    "        \n",
+    "        # Identify violations of the hit probability constraints (h > 1 or h < 0)\n",
+    "        larger_i = np.where(h_optt > 1)[0]\n",
+    "        smaller_i = np.where(h_optt < 0)[0]\n",
+    "\n",
+    "        # If no violations, the optimal solution is reached\n",
+    "        if len(smaller_i) == 0 and len(larger_i) == 0:\n",
+    "            break\n",
+    "        \n",
+    "        # Find the furthest object from the boundary (either 0 or 1)\n",
+    "        min_viol = 0\n",
+    "        min_viol_i = -1\n",
+    "        if len(smaller_i) > 0:\n",
+    "            min_viol, min_viol_i = np.min(h_optt[smaller_i]), np.argmin(h_optt[smaller_i])\n",
+    "\n",
+    "        max_viol = 0\n",
+    "        max_viol_i = -1\n",
+    "        if len(larger_i) > 0:\n",
+    "            max_viol, max_viol_i = np.max(h_optt[larger_i] - 1), np.argmax(h_optt[larger_i] - 1)\n",
+    "        \n",
+    "        # Choose the variable with the largest violation to adjust\n",
+    "        if max_viol > abs(min_viol):\n",
+    "            viol_i = max_viol_i\n",
+    "            min_viol_flag = 0\n",
+    "        else:\n",
+    "            viol_i = min_viol_i\n",
+    "            min_viol_flag = 1\n",
+    "        \n",
+    "        # Set the furthest object to the nearest boundary (0 or 1)\n",
+    "        if min_viol_flag:\n",
+    "            h_optt[viol_i] = 0\n",
+    "        else:\n",
+    "            h_optt[viol_i] = min(1, b)\n",
+    "        \n",
+    "        # Update cache size and fix the selected variable\n",
+    "        b -= h_optt[viol_i]\n",
+    "        fix_i.append(viol_i)\n",
+    "        differenc_set = np.setdiff1d(np.arange(N), fix_i)\n",
+    "        n = N - len(fix_i)\n",
+    "    \n",
+    "    return h_optt\n",
+    "\n",
+    "\n",
+    "# Example usage\n",
+    "lambda_vals = np.array(\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "5d223324-d193-416a-b2c3-1e143e981a37",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "\n",
+    "def numerical_opt(lambda_vals, B, c_f, c_delta):\n",
+    "    \"\"\"\n",
+    "    Perform numerical optimization to compute optimal hit probabilities.\n",
+    "\n",
+    "    Parameters:\n",
+    "    - lambda_vals: Array of request rates.\n",
+    "    - B: Cache size (constraint on total hit probabilities).\n",
+    "    - c_f: Cost of fetching (cache miss cost).\n",
+    "    - c_delta: Cost of caching (hit cost).\n",
+    "\n",
+    "    Returns:\n",
+    "    - x_opt: Optimal hit probabilities.\n",
+    "    \"\"\"\n",
+    "    N = len(lambda_vals)  # Number of items\n",
+    "\n",
+    "    # Initial guess: Even distribution of cache capacity\n",
+    "    x_init = np.full(N, B / N)\n",
+    "\n",
+    "    # Objective function\n",
+    "    def objective(x):\n",
+    "        return np.sum(lambda_vals * ((1 - x) * c_f + x**2 * c_delta / 2))\n",
+    "\n",
+    "    # Constraint: Sum of hit probabilities <= cache size (B)\n",
+    "    def constraint_total_hit(x):\n",
+    "        return B - np.sum(x)  # Non-negative means constraint satisfied\n",
+    "\n",
+    "    # Bounds for hit probabilities: 0 <= h_i <= 1\n",
+    "    bounds = [(0, 1) for _ in range(N)]\n",
+    "\n",
+    "    # Optimization\n",
+    "    constraints = [{'type': 'ineq', 'fun': constraint_total_hit}]  # Inequality constraint\n",
+    "    result = minimize(\n",
+    "        objective, \n",
+    "        x_init, \n",
+    "        method='SLSQP',  # Sequential Least Squares Quadratic Programming\n",
+    "        bounds=bounds, \n",
+    "        constraints=constraints, \n",
+    "        options={'disp': True}  # Set to True for optimization output\n",
+    "    )\n",
+    "\n",
+    "    # Optimal solution\n",
+    "    if result.success:\n",
+    "        return result.x  # Optimal hit probabilities\n",
+    "    else:\n",
+    "        raise ValueError(\"Optimization failed: \" + result.message)\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "073be740-dc97-454b-87e7-2f8f93c8f137",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "\n",
+    "# Example usage\n",
+    "lambda_vals = np.array([0.03, 0.04, 0.05, 0.06, 0.07, 1, 1.1, 1.2, 1.3, 1.4, 1.5])\n",
+    "B = 4.4\n",
+    "c_f = 7\n",
+    "c_delta = 1\n",
+    "\n",
+    "optimal_hit_probs = numerical_opt(lambda_vals, B, c_f, c_delta)\n",
+    "print(\"Optimal Hit Probabilities:\", optimal_hit_probs)\n",
+    "\n",
+    "\n",
+    "# Optimization\n",
+    "h_numerical = numerical_opt(lambda_vals, B, c_f, c_delta)\n",
+    "h_theoretical = theoretical_opt(lambda_vals, B, c_f, c_delta)\n",
+    "\n",
+    "# Comparison\n",
+    "hit_opt = np.vstack((h_numerical, h_theoretical)).T  # Combine for comparison\n",
+    "\n",
+    "# Objective Function Calculation\n",
+    "obj_1 = np.sum(lambda_vals * ((1 - h_numerical) * c_f + h_numerical**2 * c_delta / 2))\n",
+    "obj_2 = np.sum(lambda_vals * ((1 - h_theoretical) * c_f + h_theoretical**2 * c_delta / 2))\n",
+    "obj = [obj_1, obj_2]\n",
+    "\n",
+    "# Constraints\n",
+    "const_1 = np.sum(h_numerical) - B\n",
+    "const_2 = np.sum(h_theoretical) - B\n",
+    "constraint = [const_1, const_2]\n",
+    "\n",
+    "# Outputs\n",
+    "print(\"Hit probabilities (numerical and theoretical):\\n\", hit_opt)\n",
+    "print(\"Objective function values (numerical, theoretical):\", obj)\n",
+    "print(\"Constraint violations (numerical, theoretical):\", constraint)\n"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "graphs",
+   "language": "python",
+   "name": "graphs"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.12.7"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

01_nb_cncf_optimization/matlab/.ipynb_checkpoints/Main-checkpoint.m

@@ -0,0 +1,27 @@
+clc
+close all
+clear all
+
+
+%% Define parameters
+lambda = [0.03, 0.04,0.05,0.06,0.07,1,1.1,1.2,1.3,1.4,1.5]';  % Request rates ascendingly
+N=length(lambda);
+B = 4.4;  % Cache size
+c_delta=1; % age linear cost
+c_f=7; % fetching linear cost (caching miss cost)
+
+
+%% Optimization 
+[h_numerical ]=Numerical_opt(lambda,B,c_f,c_delta)
+[h_theo] = Theoritical_opt(lambda,B,c_f,c_delta)
+
+
+%% Comparison
+hit_opt=[h_numerical h_theo]
+obj_1=sum(lambda .* ((1-h_numerical)*c_f+h_numerical.^2*c_delta/2));
+obj_2=sum(lambda .* ((1-h_theo)*c_f+h_theo.^2*c_delta/2));
+obj=[obj_1 obj_2]
+const_1=sum(h_numerical)-B;
+const_2=sum(h_theo)-B;
+constraint=[const_1 const_2]
+

01_nb_cncf_optimization/matlab/.ipynb_checkpoints/Numerical_opt-checkpoint.m

@@ -0,0 +1,17 @@
+function [x_opt] = Numerical_opt(lambda,B,c_f,c_delta)
+% Numerical optimization
+% x_opt are the optimal hit probabilities
+ 
+N = length(lambda);  % Number of variables
+h_init = ones(N, 1) * B / N;  % Initial guess for h_i, evenly distributed
+x_init=h_init;
+objective1 = @(x) sum(lambda .*( (1-x)*c_f+x.^2*c_delta/2));  % Objective function
+constraint1 = @(x) sum(x) - B; % cache size constraint
+% constraint on hit prob. 0<=h_i<=1
+nonlcon1 = @(x) deal(constraint1(x), []);  % 
+lb = zeros(N, 1);  % Lower bounds
+ub = ones(N, 1);     % Upper bounds
+options = optimoptions('fmincon', 'Display', 'iter', 'Algorithm', 'sqp');
+[x_opt, fval_h] = fmincon(objective1, x_init, [], [], [], [], lb, ub, nonlcon1, options);
+end
+

01_nb_cncf_optimization/matlab/.ipynb_checkpoints/Theoritical_opt-Copy1-checkpoint.m

@@ -0,0 +1,171 @@
+function [h_optt] = Theoritical_opt(lambda,B,c_f,c_delta)
+%% Theoritical optimization 
+
+
+%% Iterative identification of active constraints
+N=length(lambda)
+flag=1;
+h_optt=zeros(N,1); %optimal hit prob
+differenc_set=1:N; % the set of variables to optimize 
+fix_i=[]; % set of variables that reached optimality and are excluded from the optimization
+n=N; 
+b=B;
+
+%%
+while flag
+    if(n==0) 
+        if(b>0) % if there is left over cache size and mu is not zero (the loop would break), redistribute it among the zero hit probability
+            differenc_set=find(h_optt==0)';
+            fix_i=setdiff(1:N,differenc_set)';
+            n=length(differenc_set);
+            continue;     
+        else           
+             h_optt(differenc_set)=0;
+             break;
+        end
+    end
+    % Optimal Lagrangian mult. and hit prob. calculated theoritically for the set of variables in  differenc_set
+    mu=max(0,(n*c_f-b*c_delta)/ sum(1./lambda(differenc_set)));  %optimal lagrangian mult.
+    h_optt(differenc_set)=(c_f-mu./lambda(differenc_set))/c_delta %optimal hit prob
+    
+    % mu has to be >=0
+    if(mu<0)
+        b=(n*c_f/c_delta); % this sets mu to zero in the next iteration
+        continue;
+    end
+    
+    % check the violation of the hit_prob const  
+    larger_i=find(h_optt>1); % h>1
+    smaller_i=find(h_optt<0); % h<0
+
+%     smaller=h(find(differenc_set<0))-0;
+    % no violation means optimal solution is reached for all objects
+    if(length(smaller_i)==0 && length(larger_i)==0)
+%         flag=0;
+  
+        break;
+    end
+    
+    % find the furthest object from the 0 boundary
+    min_viol=0;
+    min_viol_i=-1;
+    if(length(smaller_i)>0)
+        [min_viol, min_viol_i]=min(h_optt);
+    end
+     % find the furthest object from the 1 boundary
+    max_viol=0; 
+    max_viol_i=-1;
+    if(length(larger_i)>0)
+       larger=h_optt-1;
+       [max_viol ,max_viol_i]=max(h_optt-1);
+
+    end
+    
+     % compare both furthest objects from both boundaries
+    viol_i=min_viol_i;
+    min_viol_flag=1; % True if the furthest one is from the left
+    if(max_viol>abs(min_viol))
+        viol_i= max_viol_i;
+        min_viol_flag=0;
+
+    end
+     % set the furthest object to the nearest boundary
+     if(min_viol_flag)
+        h_optt(viol_i)=0;
+    else
+        h_optt(viol_i)=min(1,b);
+     end
+    
+    %calculate the new parameters after removing the furthest object from
+    %the decision variables
+    B_new=b-(h_optt(viol_i));
+    b=B_new;
+    fix_i=[fix_i' viol_i']';
+    differenc_set=setdiff(1:N,fix_i)   ; 
+    n=N-length(fix_i);
+    
+    
+    %     % Identify the most violating object from the right side  h>1
+%     if(length(larger_i)>0)
+%        larger=h_optt-1;
+%        [max_viol ,max_viol_i]=max(h_optt-1); % maximum violating object
+%        h_optt(max_viol_i)=min(1,b); %project to the feasible range
+%        b=max(b-1,0); % update the cache size 
+%        fix_i=[fix_i' max_viol_i']'; %exclude i from the set of decision variables
+%        differenc_set=setdiff(1:N,fix_i); % obtain the set of decision variables 
+%        n=N-length(fix_i); % update the number of decision variables 
+%        continue;
+%     end
+%     
+%     if(length(smaller_i)>0)
+%         [min_viol, min_viol_i]=min(h_optt);
+%         h_optt(min_viol_i)=0; 
+%         fix_i=[fix_i' min_viol_i']';
+%         differenc_set=setdiff(1:N,fix_i)   ; 
+%         n=N-length(fix_i);
+%     end
+%     
+% end
+
+end
+
+%% Identfying the active constraints collectively
+% flag=1;
+% h_optt=zeros(N,1);
+% differenc_set=1:N;
+% fix_i=[];
+% n=N;
+% b=B;
+% while flag
+%     mu=(n*c_f-b*c_delta)/ sum(1./lambda(differenc_set));
+%     h_optt(differenc_set)=(c_f-mu./lambda(differenc_set))/c_delta
+% 
+%     larger_i=find(h_optt>1);    
+% %     larger=h_optt(larger_i)-1;
+%     smaller_i=find(h_optt<0);
+% %     smaller=h(find(differenc_set<0))-0;
+%     mult=solve_multipliers(lambda,B,c_f,c_delta,larger_i)
+% 
+%     if(length(larger_i)+length(smaller_i)==0)
+%         flag=0;
+%         break;
+%     end
+%     if(length(smaller_i)>0)
+%         h_optt(smaller_i)=0;
+%         fix_i=[fix_i' smaller_i' ]';
+%         differenc_set=setdiff(1:N,fix_i)    
+%         n=N-length(fix_i);
+%         continue
+%     end
+% %     h_optt(smaller_i)=0;
+%     if(length(larger_i)>b)
+%         [~,index]=maxk(h_optt,b)
+%         h_optt(index)=1;
+%         B_new=b-sum(h_optt(index));
+%         fix_i=[fix_i' smaller_i' index']';
+%     else
+%          h_optt(larger_i)=1;
+%          B_new=b-sum(h_optt(larger_i));
+%          fix_i=[fix_i' smaller_i' larger_i']';
+%     end
+% %     mult=solve_multipliers(lambda,B,c_f,c_delta,larger_i)
+% 
+%     b=B_new;
+%    
+%     differenc_set=setdiff(1:N,fix_i)    
+%     n=N-length(fix_i);
+% end
+% % h_optt=zeros(N,1);
+% % h_optt(end-B+1:end)=1;
+% optimal=[h_opt h_optt]
+% obj_1=sum(lambda .* ((1-h_opt)*c_f+h_opt.^2*c_delta/2));
+% obj_2=sum(lambda .* ((1-h_optt)*c_f+h_optt.^2*c_delta/2));
+% objective=[obj_1 obj_2]
+% const_1=sum(h_opt)-B;
+% const_2=sum(h_optt)-B;
+% constraint=[const_1 const_2]
+% 
+
+
+end
+

01_nb_cncf_optimization/matlab/.ipynb_checkpoints/Theoritical_opt-checkpoint.m

@@ -0,0 +1,171 @@
+function [h_optt] = Theoritical_opt(lambda,B,c_f,c_delta)
+%% Theoritical optimization 
+
+
+%% Iterative identification of active constraints
+N=length(lambda)
+flag=1;
+h_optt=zeros(N,1); %optimal hit prob
+differenc_set=1:N; % the set of variables to optimize 
+fix_i=[]; % set of variables that reached optimality and are excluded from the optimization
+n=N; 
+b=B;
+
+%%
+while flag
+    if(n==0) 
+        if(b>0) % if there is left over cache size and mu is not zero (the loop would break), redistribute it among the zero hit probability
+            differenc_set=find(h_optt==0)';
+            fix_i=setdiff(1:N,differenc_set)';
+            n=length(differenc_set);
+            continue;     
+        else           
+             h_optt(differenc_set)=0;
+             break;
+        end
+    end
+    % Optimal Lagrangian mult. and hit prob. calculated theoritically for the set of variables in  differenc_set
+    mu=max(0,(n*c_f-b*c_delta)/ sum(1./lambda(differenc_set)));  %optimal lagrangian mult.
+    h_optt(differenc_set)=(c_f-mu./lambda(differenc_set))/c_delta %optimal hit prob
+    
+    % mu has to be >=0
+    if(mu<0)
+        b=(n*c_f/c_delta); % this sets mu to zero in the next iteration
+        continue;
+    end
+    
+    % check the violation of the hit_prob const  
+    larger_i=find(h_optt>1); % h>1
+    smaller_i=find(h_optt<0); % h<0
+
+%     smaller=h(find(differenc_set<0))-0;
+    % no violation means optimal solution is reached for all objects
+    if(length(smaller_i)==0 && length(larger_i)==0)
+%         flag=0;
+  
+        break;
+    end
+    
+    % find the furthest object from the 0 boundary
+    min_viol=0;
+    min_viol_i=-1;
+    if(length(smaller_i)>0)
+        [min_viol, min_viol_i]=min(h_optt);
+    end
+     % find the furthest object from the 1 boundary
+    max_viol=0; 
+    max_viol_i=-1;
+    if(length(larger_i)>0)
+       larger=h_optt-1;
+       [max_viol ,max_viol_i]=max(h_optt-1);
+
+    end
+    
+     % compare both furthest objects from both boundaries
+    viol_i=min_viol_i;
+    min_viol_flag=1; % True if the furthest one is from the left
+    if(max_viol>abs(min_viol))
+        viol_i= max_viol_i;
+        min_viol_flag=0;
+
+    end
+     % set the furthest object to the nearest boundary
+     if(min_viol_flag)
+        h_optt(viol_i)=0;
+    else
+        h_optt(viol_i)=min(1,b);
+     end
+    
+    %calculate the new parameters after removing the furthest object from
+    %the decision variables
+    B_new=b-(h_optt(viol_i));
+    b=B_new;
+    fix_i=[fix_i' viol_i']';
+    differenc_set=setdiff(1:N,fix_i)   ; 
+    n=N-length(fix_i);
+    
+    
+    %     % Identify the most violating object from the right side  h>1
+%     if(length(larger_i)>0)
+%        larger=h_optt-1;
+%        [max_viol ,max_viol_i]=max(h_optt-1); % maximum violating object
+%        h_optt(max_viol_i)=min(1,b); %project to the feasible range
+%        b=max(b-1,0); % update the cache size 
+%        fix_i=[fix_i' max_viol_i']'; %exclude i from the set of decision variables
+%        differenc_set=setdiff(1:N,fix_i); % obtain the set of decision variables 
+%        n=N-length(fix_i); % update the number of decision variables 
+%        continue;
+%     end
+%     
+%     if(length(smaller_i)>0)
+%         [min_viol, min_viol_i]=min(h_optt);
+%         h_optt(min_viol_i)=0; 
+%         fix_i=[fix_i' min_viol_i']';
+%         differenc_set=setdiff(1:N,fix_i)   ; 
+%         n=N-length(fix_i);
+%     end
+%     
+% end
+
+end
+
+%% Identfying the active constraints collectively
+% flag=1;
+% h_optt=zeros(N,1);
+% differenc_set=1:N;
+% fix_i=[];
+% n=N;
+% b=B;
+% while flag
+%     mu=(n*c_f-b*c_delta)/ sum(1./lambda(differenc_set));
+%     h_optt(differenc_set)=(c_f-mu./lambda(differenc_set))/c_delta
+% 
+%     larger_i=find(h_optt>1);    
+% %     larger=h_optt(larger_i)-1;
+%     smaller_i=find(h_optt<0);
+% %     smaller=h(find(differenc_set<0))-0;
+%     mult=solve_multipliers(lambda,B,c_f,c_delta,larger_i)
+% 
+%     if(length(larger_i)+length(smaller_i)==0)
+%         flag=0;
+%         break;
+%     end
+%     if(length(smaller_i)>0)
+%         h_optt(smaller_i)=0;
+%         fix_i=[fix_i' smaller_i' ]';
+%         differenc_set=setdiff(1:N,fix_i)    
+%         n=N-length(fix_i);
+%         continue
+%     end
+% %     h_optt(smaller_i)=0;
+%     if(length(larger_i)>b)
+%         [~,index]=maxk(h_optt,b)
+%         h_optt(index)=1;
+%         B_new=b-sum(h_optt(index));
+%         fix_i=[fix_i' smaller_i' index']';
+%     else
+%          h_optt(larger_i)=1;
+%          B_new=b-sum(h_optt(larger_i));
+%          fix_i=[fix_i' smaller_i' larger_i']';
+%     end
+% %     mult=solve_multipliers(lambda,B,c_f,c_delta,larger_i)
+% 
+%     b=B_new;
+%    
+%     differenc_set=setdiff(1:N,fix_i)    
+%     n=N-length(fix_i);
+% end
+% % h_optt=zeros(N,1);
+% % h_optt(end-B+1:end)=1;
+% optimal=[h_opt h_optt]
+% obj_1=sum(lambda .* ((1-h_opt)*c_f+h_opt.^2*c_delta/2));
+% obj_2=sum(lambda .* ((1-h_optt)*c_f+h_optt.^2*c_delta/2));
+% objective=[obj_1 obj_2]
+% const_1=sum(h_opt)-B;
+% const_2=sum(h_optt)-B;
+% constraint=[const_1 const_2]
+% 
+
+
+end
+

01_nb_cncf_optimization/matlab/Main.ipynb

@@ -0,0 +1,248 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "id": "d0996120-bb17-4476-b912-ce155100b2cb",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "from scipy.optimize import minimize\n",
+    "import numpy as np\n",
+    "\n",
+    "# Define Parameters\n",
+    "lambda_vals = np.array([0.03, 0.04, 0.05, 0.06, 0.07, 1, 1.1, 1.2, 1.3, 1.4, 1.5])  # Request rates ascendingly\n",
+    "N = len(lambda_vals)\n",
+    "B = 4.4  # Cache size\n",
+    "c_delta = 1  # Age linear cost\n",
+    "c_f = 7  # Fetching linear cost (caching miss cost)\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "id": "fcf0c13c-5b2c-457e-9aa6-8d349fcf13fa",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "\n",
+    "def theoretical_opt(lambda_vals, B, c_f, c_delta):\n",
+    "    \"\"\"\n",
+    "    Perform theoretical optimization to compute optimal hit probabilities.\n",
+    "    \n",
+    "    Parameters:\n",
+    "    - lambda_vals: Array of request rates.\n",
+    "    - B: Cache size (constraint on total hit probabilities).\n",
+    "    - c_f: Cost of fetching (cache miss cost).\n",
+    "    - c_delta: Cost of caching (hit cost).\n",
+    "\n",
+    "    Returns:\n",
+    "    - h_optt: Optimal hit probabilities.\n",
+    "    \"\"\"\n",
+    "    N = len(lambda_vals)\n",
+    "    h_optt = np.zeros(N)  # Initialize optimal hit probabilities\n",
+    "    differenc_set = np.arange(N)  # Set of variables to optimize\n",
+    "    fix_i = []  # Set of fixed variables (those already optimized)\n",
+    "    n = N\n",
+    "    b = B\n",
+    "    flag = True\n",
+    "\n",
+    "    while flag:\n",
+    "        if n == 0:  # If no variables left to optimize\n",
+    "            if b > 0:  # If there is leftover cache size, redistribute it\n",
+    "                differenc_set = np.where(h_optt == 0)[0]  # Find zero hit probability variables\n",
+    "                fix_i = np.setdiff1d(np.arange(N), differenc_set)\n",
+    "                n = len(differenc_set)\n",
+    "                continue\n",
+    "            else:  # No variables to optimize, finalize\n",
+    "                h_optt[differenc_set] = 0\n",
+    "                break\n",
+    "        \n",
+    "        # Calculate the optimal Lagrange multiplier (mu) and hit probabilities for the set of variables\n",
+    "        mu = max(0, (n * c_f - b * c_delta) / np.sum(1 / lambda_vals[differenc_set]))\n",
+    "        h_optt[differenc_set] = (c_f - mu / lambda_vals[differenc_set]) / c_delta\n",
+    "        \n",
+    "        # If mu < 0, adjust the cache size to set mu to zero in the next iteration\n",
+    "        if mu < 0:\n",
+    "            b = (n * c_f / c_delta)\n",
+    "            continue\n",
+    "        \n",
+    "        # Identify violations of the hit probability constraints (h > 1 or h < 0)\n",
+    "        larger_i = np.where(h_optt > 1)[0]\n",
+    "        smaller_i = np.where(h_optt < 0)[0]\n",
+    "\n",
+    "        # If no violations, the optimal solution is reached\n",
+    "        if len(smaller_i) == 0 and len(larger_i) == 0:\n",
+    "            break\n",
+    "        \n",
+    "        # Find the furthest object from the boundary (either 0 or 1)\n",
+    "        min_viol = 0\n",
+    "        min_viol_i = -1\n",
+    "        if len(smaller_i) > 0:\n",
+    "            min_viol, min_viol_i = np.min(h_optt[smaller_i]), np.argmin(h_optt[smaller_i])\n",
+    "\n",
+    "        max_viol = 0\n",
+    "        max_viol_i = -1\n",
+    "        if len(larger_i) > 0:\n",
+    "            max_viol, max_viol_i = np.max(h_optt[larger_i] - 1), np.argmax(h_optt[larger_i] - 1)\n",
+    "        \n",
+    "        # Choose the variable with the largest violation to adjust\n",
+    "        if max_viol > abs(min_viol):\n",
+    "            viol_i = max_viol_i\n",
+    "            min_viol_flag = 0\n",
+    "        else:\n",
+    "            viol_i = min_viol_i\n",
+    "            min_viol_flag = 1\n",
+    "        \n",
+    "        # Set the furthest object to the nearest boundary (0 or 1)\n",
+    "        if min_viol_flag:\n",
+    "            h_optt[viol_i] = 0\n",
+    "        else:\n",
+    "            h_optt[viol_i] = min(1, b)\n",
+    "        \n",
+    "        # Update cache size and fix the selected variable\n",
+    "        b -= h_optt[viol_i]\n",
+    "        fix_i.append(viol_i)\n",
+    "        differenc_set = np.setdiff1d(np.arange(N), fix_i)\n",
+    "        n = N - len(fix_i)\n",
+    "    \n",
+    "    return h_optt"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "id": "5d223324-d193-416a-b2c3-1e143e981a37",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "\n",
+    "def numerical_opt(lambda_vals, B, c_f, c_delta):\n",
+    "    \"\"\"\n",
+    "    Perform numerical optimization to compute optimal hit probabilities.\n",
+    "\n",
+    "    Parameters:\n",
+    "    - lambda_vals: Array of request rates.\n",
+    "    - B: Cache size (constraint on total hit probabilities).\n",
+    "    - c_f: Cost of fetching (cache miss cost).\n",
+    "    - c_delta: Cost of caching (hit cost).\n",
+    "\n",
+    "    Returns:\n",
+    "    - x_opt: Optimal hit probabilities.\n",
+    "    \"\"\"\n",
+    "    N = len(lambda_vals)  # Number of items\n",
+    "\n",
+    "    # Initial guess: Even distribution of cache capacity\n",
+    "    x_init = np.full(N, B / N)\n",
+    "\n",
+    "    # Objective function\n",
+    "    def objective(x):\n",
+    "        return np.sum(lambda_vals * ((1 - x) * c_f + x**2 * c_delta / 2))\n",
+    "\n",
+    "    # Constraint: Sum of hit probabilities <= cache size (B)\n",
+    "    def constraint_total_hit(x):\n",
+    "        return B - np.sum(x)  # Non-negative means constraint satisfied\n",
+    "\n",
+    "    # Bounds for hit probabilities: 0 <= h_i <= 1\n",
+    "    bounds = [(0, 1) for _ in range(N)]\n",
+    "\n",
+    "    # Optimization\n",
+    "    constraints = [{'type': 'ineq', 'fun': constraint_total_hit}]  # Inequality constraint\n",
+    "    result = minimize(\n",
+    "        objective, \n",
+    "        x_init, \n",
+    "        method='SLSQP',  # Sequential Least Squares Quadratic Programming\n",
+    "        bounds=bounds, \n",
+    "        constraints=constraints, \n",
+    "        options={'disp': True}  # Set to True for optimization output\n",
+    "    )\n",
+    "\n",
+    "    # Optimal solution\n",
+    "    if result.success:\n",
+    "        return result.x  # Optimal hit probabilities\n",
+    "    else:\n",
+    "        raise ValueError(\"Optimization failed: \" + result.message)\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "id": "073be740-dc97-454b-87e7-2f8f93c8f137",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Optimization terminated successfully    (Exit mode 0)\n",
+      "            Current function value: 16.15721739129548\n",
+      "            Iterations: 4\n",
+      "            Function evaluations: 48\n",
+      "            Gradient evaluations: 4\n",
+      "Hit probabilities (numerical and theoretical):\n",
+      " [[0.00000000e+00 0.00000000e+00]\n",
+      " [0.00000000e+00 0.00000000e+00]\n",
+      " [0.00000000e+00 0.00000000e+00]\n",
+      " [0.00000000e+00 0.00000000e+00]\n",
+      " [0.00000000e+00 0.00000000e+00]\n",
+      " [4.58923826e-13 1.00000000e+00]\n",
+      " [4.26087031e-01 0.00000000e+00]\n",
+      " [9.73912969e-01 4.00000000e-01]\n",
+      " [1.00000000e+00 0.00000000e+00]\n",
+      " [1.00000000e+00 0.00000000e+00]\n",
+      " [1.00000000e+00 0.00000000e+00]]\n",
+      "Objective function values (numerical, theoretical): [16.15721739129548, 44.486]\n",
+      "Constraint violations (numerical, theoretical): [1.241673430740775e-12, -3.0]\n"
+     ]
+    }
+   ],
+   "source": [
+    "\n",
+    "# Optimization\n",
+    "h_numerical = numerical_opt(lambda_vals, B, c_f, c_delta)\n",
+    "h_theoretical = theoretical_opt(lambda_vals, B, c_f, c_delta)\n",
+    "\n",
+    "# Comparison of Hit Probabilities\n",
+    "hit_opt = np.vstack((h_numerical, h_theoretical)).T  # Combine numerical and theoretical hit probabilities\n",
+    "\n",
+    "# Objective Function Calculation\n",
+    "obj_1 = np.sum(lambda_vals * ((1 - h_numerical) * c_f + h_numerical**2 * c_delta / 2))\n",
+    "obj_2 = np.sum(lambda_vals * ((1 - h_theoretical) * c_f + h_theoretical**2 * c_delta / 2))\n",
+    "obj = [obj_1, obj_2]  # Store objective function values for both methods\n",
+    "\n",
+    "# Constraints\n",
+    "const_1 = np.sum(h_numerical) - B\n",
+    "const_2 = np.sum(h_theoretical) - B\n",
+    "constraint = [const_1, const_2]  # Check if the cache size constraint is satisfied\n",
+    "\n",
+    "# Outputs\n",
+    "print(\"Hit probabilities (numerical and theoretical):\\n\", hit_opt)\n",
+    "print(\"Objective function values (numerical, theoretical):\", obj)\n",
+    "print(\"Constraint violations (numerical, theoretical):\", constraint)"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "graphs",
+   "language": "python",
+   "name": "graphs"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.12.7"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

01_nb_cncf_optimization/matlab/Main.m

@@ -0,0 +1,27 @@
+clc
+close all
+clear all
+
+
+%% Define parameters
+lambda = [0.03, 0.04,0.05,0.06,0.07,1,1.1,1.2,1.3,1.4,1.5]';  % Request rates ascendingly
+N=length(lambda);
+B = 4.4;  % Cache size
+c_delta=1; % age linear cost
+c_f=7; % fetching linear cost (caching miss cost)
+
+
+%% Optimization 
+[h_numerical ]=Numerical_opt(lambda,B,c_f,c_delta)
+[h_theo] = Theoritical_opt(lambda,B,c_f,c_delta)
+
+
+%% Comparison
+hit_opt=[h_numerical h_theo]
+obj_1=sum(lambda .* ((1-h_numerical)*c_f+h_numerical.^2*c_delta/2));
+obj_2=sum(lambda .* ((1-h_theo)*c_f+h_theo.^2*c_delta/2));
+obj=[obj_1 obj_2]
+const_1=sum(h_numerical)-B;
+const_2=sum(h_theo)-B;
+constraint=[const_1 const_2]
+

01_nb_cncf_optimization/matlab/Numerical_opt.m

@@ -0,0 +1,17 @@
+function [x_opt] = Numerical_opt(lambda,B,c_f,c_delta)
+% Numerical optimization
+% x_opt are the optimal hit probabilities
+ 
+N = length(lambda);  % Number of variables
+h_init = ones(N, 1) * B / N;  % Initial guess for h_i, evenly distributed
+x_init=h_init;
+objective1 = @(x) sum(lambda .*( (1-x)*c_f+x.^2*c_delta/2));  % Objective function
+constraint1 = @(x) sum(x) - B; % cache size constraint
+% constraint on hit prob. 0<=h_i<=1
+nonlcon1 = @(x) deal(constraint1(x), []);  % 
+lb = zeros(N, 1);  % Lower bounds
+ub = ones(N, 1);     % Upper bounds
+options = optimoptions('fmincon', 'Display', 'iter', 'Algorithm', 'sqp');
+[x_opt, fval_h] = fmincon(objective1, x_init, [], [], [], [], lb, ub, nonlcon1, options);
+end
+

01_nb_cncf_optimization/matlab/Theoritical_opt-Copy1.m

@@ -0,0 +1,90 @@
+function [h_optt] = Theoritical_opt(lambda,B,c_f,c_delta)
+%% Theoritical optimization 
+
+
+%% Iterative identification of active constraints
+N=length(lambda)
+flag=1;
+h_optt=zeros(N,1); %optimal hit prob
+differenc_set=1:N; % the set of variables to optimize 
+fix_i=[]; % set of variables that reached optimality and are excluded from the optimization
+n=N; 
+b=B;
+
+%%
+while flag
+    if(n==0) 
+        if(b>0) % if there is left over cache size and mu is not zero (the loop would break), redistribute it among the zero hit probability
+            differenc_set=find(h_optt==0)';
+            fix_i=setdiff(1:N,differenc_set)';
+            n=length(differenc_set);
+            continue;     
+        else           
+             h_optt(differenc_set)=0;
+             break;
+        end
+    end
+    % Optimal Lagrangian mult. and hit prob. calculated theoritically for the set of variables in  differenc_set
+    mu=max(0,(n*c_f-b*c_delta)/ sum(1./lambda(differenc_set)));  %optimal lagrangian mult.
+    h_optt(differenc_set)=(c_f-mu./lambda(differenc_set))/c_delta %optimal hit prob
+    
+    % mu has to be >=0
+    if(mu<0)
+        b=(n*c_f/c_delta); % this sets mu to zero in the next iteration
+        continue;
+    end
+    
+    % check the violation of the hit_prob const  
+    larger_i=find(h_optt>1); % h>1
+    smaller_i=find(h_optt<0); % h<0
+
+    if(length(smaller_i)==0 && length(larger_i)==0)
+
+  
+        break;
+    end
+    
+    % find the furthest object from the 0 boundary
+    min_viol=0;
+    min_viol_i=-1;
+    if(length(smaller_i)>0)
+        [min_viol, min_viol_i]=min(h_optt);
+    end
+     % find the furthest object from the 1 boundary
+    max_viol=0; 
+    max_viol_i=-1;
+    if(length(larger_i)>0)
+       larger=h_optt-1;
+       [max_viol ,max_viol_i]=max(h_optt-1);
+
+    end
+    
+     % compare both furthest objects from both boundaries
+    viol_i=min_viol_i;
+    min_viol_flag=1; % True if the furthest one is from the left
+    if(max_viol>abs(min_viol))
+        viol_i= max_viol_i;
+        min_viol_flag=0;
+
+    end
+     % set the furthest object to the nearest boundary
+     if(min_viol_flag)
+        h_optt(viol_i)=0;
+    else
+        h_optt(viol_i)=min(1,b);
+     end
+    
+    %calculate the new parameters after removing the furthest object from
+    %the decision variables
+    B_new=b-(h_optt(viol_i));
+    b=B_new;
+    fix_i=[fix_i' viol_i']';
+    differenc_set=setdiff(1:N,fix_i)   ; 
+    n=N-length(fix_i);
+    
+
+
+end
+
+end
+

01_nb_cncf_optimization/matlab/Theoritical_opt.m

@@ -0,0 +1,171 @@
+function [h_optt] = Theoritical_opt(lambda,B,c_f,c_delta)
+%% Theoritical optimization 
+
+
+%% Iterative identification of active constraints
+N=length(lambda)
+flag=1;
+h_optt=zeros(N,1); %optimal hit prob
+differenc_set=1:N; % the set of variables to optimize 
+fix_i=[]; % set of variables that reached optimality and are excluded from the optimization
+n=N; 
+b=B;
+
+%%
+while flag
+    if(n==0) 
+        if(b>0) % if there is left over cache size and mu is not zero (the loop would break), redistribute it among the zero hit probability
+            differenc_set=find(h_optt==0)';
+            fix_i=setdiff(1:N,differenc_set)';
+            n=length(differenc_set);
+            continue;     
+        else           
+             h_optt(differenc_set)=0;
+             break;
+        end
+    end
+    % Optimal Lagrangian mult. and hit prob. calculated theoritically for the set of variables in  differenc_set
+    mu=max(0,(n*c_f-b*c_delta)/ sum(1./lambda(differenc_set)));  %optimal lagrangian mult.
+    h_optt(differenc_set)=(c_f-mu./lambda(differenc_set))/c_delta %optimal hit prob
+    
+    % mu has to be >=0
+    if(mu<0)
+        b=(n*c_f/c_delta); % this sets mu to zero in the next iteration
+        continue;
+    end
+    
+    % check the violation of the hit_prob const  
+    larger_i=find(h_optt>1); % h>1
+    smaller_i=find(h_optt<0); % h<0
+
+%     smaller=h(find(differenc_set<0))-0;
+    % no violation means optimal solution is reached for all objects
+    if(length(smaller_i)==0 && length(larger_i)==0)
+%         flag=0;
+  
+        break;
+    end
+    
+    % find the furthest object from the 0 boundary
+    min_viol=0;
+    min_viol_i=-1;
+    if(length(smaller_i)>0)
+        [min_viol, min_viol_i]=min(h_optt);
+    end
+     % find the furthest object from the 1 boundary
+    max_viol=0; 
+    max_viol_i=-1;
+    if(length(larger_i)>0)
+       larger=h_optt-1;
+       [max_viol ,max_viol_i]=max(h_optt-1);
+
+    end
+    
+     % compare both furthest objects from both boundaries
+    viol_i=min_viol_i;
+    min_viol_flag=1; % True if the furthest one is from the left
+    if(max_viol>abs(min_viol))
+        viol_i= max_viol_i;
+        min_viol_flag=0;
+
+    end
+     % set the furthest object to the nearest boundary
+     if(min_viol_flag)
+        h_optt(viol_i)=0;
+    else
+        h_optt(viol_i)=min(1,b);
+     end
+    
+    %calculate the new parameters after removing the furthest object from
+    %the decision variables
+    B_new=b-(h_optt(viol_i));
+    b=B_new;
+    fix_i=[fix_i' viol_i']';
+    differenc_set=setdiff(1:N,fix_i)   ; 
+    n=N-length(fix_i);
+    
+    
+    %     % Identify the most violating object from the right side  h>1
+%     if(length(larger_i)>0)
+%        larger=h_optt-1;
+%        [max_viol ,max_viol_i]=max(h_optt-1); % maximum violating object
+%        h_optt(max_viol_i)=min(1,b); %project to the feasible range
+%        b=max(b-1,0); % update the cache size 
+%        fix_i=[fix_i' max_viol_i']'; %exclude i from the set of decision variables
+%        differenc_set=setdiff(1:N,fix_i); % obtain the set of decision variables 
+%        n=N-length(fix_i); % update the number of decision variables 
+%        continue;
+%     end
+%     
+%     if(length(smaller_i)>0)
+%         [min_viol, min_viol_i]=min(h_optt);
+%         h_optt(min_viol_i)=0; 
+%         fix_i=[fix_i' min_viol_i']';
+%         differenc_set=setdiff(1:N,fix_i)   ; 
+%         n=N-length(fix_i);
+%     end
+%     
+% end
+
+end
+
+%% Identfying the active constraints collectively
+% flag=1;
+% h_optt=zeros(N,1);
+% differenc_set=1:N;
+% fix_i=[];
+% n=N;
+% b=B;
+% while flag
+%     mu=(n*c_f-b*c_delta)/ sum(1./lambda(differenc_set));
+%     h_optt(differenc_set)=(c_f-mu./lambda(differenc_set))/c_delta
+% 
+%     larger_i=find(h_optt>1);    
+% %     larger=h_optt(larger_i)-1;
+%     smaller_i=find(h_optt<0);
+% %     smaller=h(find(differenc_set<0))-0;
+%     mult=solve_multipliers(lambda,B,c_f,c_delta,larger_i)
+% 
+%     if(length(larger_i)+length(smaller_i)==0)
+%         flag=0;
+%         break;
+%     end
+%     if(length(smaller_i)>0)
+%         h_optt(smaller_i)=0;
+%         fix_i=[fix_i' smaller_i' ]';
+%         differenc_set=setdiff(1:N,fix_i)    
+%         n=N-length(fix_i);
+%         continue
+%     end
+% %     h_optt(smaller_i)=0;
+%     if(length(larger_i)>b)
+%         [~,index]=maxk(h_optt,b)
+%         h_optt(index)=1;
+%         B_new=b-sum(h_optt(index));
+%         fix_i=[fix_i' smaller_i' index']';
+%     else
+%          h_optt(larger_i)=1;
+%          B_new=b-sum(h_optt(larger_i));
+%          fix_i=[fix_i' smaller_i' larger_i']';
+%     end
+% %     mult=solve_multipliers(lambda,B,c_f,c_delta,larger_i)
+% 
+%     b=B_new;
+%    
+%     differenc_set=setdiff(1:N,fix_i)    
+%     n=N-length(fix_i);
+% end
+% % h_optt=zeros(N,1);
+% % h_optt(end-B+1:end)=1;
+% optimal=[h_opt h_optt]
+% obj_1=sum(lambda .* ((1-h_opt)*c_f+h_opt.^2*c_delta/2));
+% obj_2=sum(lambda .* ((1-h_optt)*c_f+h_optt.^2*c_delta/2));
+% objective=[obj_1 obj_2]
+% const_1=sum(h_opt)-B;
+% const_2=sum(h_optt)-B;
+% constraint=[const_1 const_2]
+% 
+
+
+end
+