Optimization of RC Beam under Flexure using Harmony Search (ACI 318-19)

In this blog, we present an example of RC beam section optimization. The optimization problem is formulated as following:

Design Vars: b, d, rho

Min cost (sum of concrete and rebar)

s.t.

phi*Mn >= Mu (given)

rho_min (0.003) <= rho <= rho_max (0.025)

bound on b,d

% Example: RC beam (rectangular, singly reinforced), flexural design
% Variables x = [b, d, rho]
%   b   : width (m)
%   d   : effective depth (m)
%   rho : steel ratio As/(b*d)
%
% Objective: minimize cost = Cc * b*h + Cs * As
%   with h = d + cover (approx), As = rho*b*d
%
% Constraints:
%   1) phi*Mn >= Mu, where:
%       a   = As*fy / (0.85*fc*b)
%       Mn  = As*fy*(d - a/2)
%       phi = 0.9 (tension-controlled, simplified)
%   2) rho_min <= rho <= rho_max
%   3) bounds on b, d
% Optional: add serviceability, max bar spacing, etc.

clear; clc;

% Known data (SI)
fc   = 30e6;          % Pa
fy   = 500e6;         % Pa
Mu   = 250e3;         % N·m (example factored moment)
phi  = 0.9;
gamma_c = 24e3;       % N/m^3 (not used directly)
cover = 0.06;         % m (approx; for cost)
Cc   = 2500;          % THB/m^3 concrete unit cost (example)
Cs   = 40000;         % THB/ton steel (example) -> ~ THB/kg; convert via density
rho_min = 0.003;      % example min steel ratio
rho_max = 0.025;      % example max steel ratio

% Objective
prob.obj = @(x) obj_cost(x, Cc, Cs, cover);   % see subfunction

% Decision variable bounds [b, d, rho]
prob.lb = [0.25, 0.40, rho_min];
prob.ub = [0.80, 0.90, rho_max];

% No integers here (could add bar-count as integer later)
prob.isInt = [false, false, false];

% No linear constraints in this simple example
prob.A = []; prob.b = []; prob.Aeq = []; prob.beq = [];

% Nonlinear constraints (explicit functions)
prob.nl_ineq = @(x) nl_ineq_rc(x, fc, fy, Mu, phi);  % g(x)<=0
prob.nl_eq   = [];                                    % none

% Optional repair (e.g., enforce practical d >= b, etc.)
prob.repair = @(x) repair_rc(x);

% HS options
opts = struct('HMS',30,'HMCR',0.9,'PAR',0.3,'bw',[0.05 0.05 0.002], ...
              'bwDamp',0.995,'maxIter',3000,'seed',2);

out = hs_opt_general(prob, opts);

disp('Best design (b, d, rho):');
disp(out.best.x);
fprintf('Best cost: %.2f\n', out.best.f);
fprintf('Feasible: %d; Violation %.3e\n', out.best.feasible, out.best.vio);

figure; 
yyaxis left; plot(out.hist.bestf,'-'); ylabel('best cost');
yyaxis right; plot(out.hist.bestvio,'--'); ylabel('best vio'); grid on;
xlabel('iteration'); title('HS on RC Beam Flexural Design');

Objective function

% -------------------- subfunctions --------------------
function f = obj_cost(x, Cc, Cs, cover)
    b = x(1); d = x(2); rho = x(3);
    h = d + cover;
    As = rho*b*d;                % m^2
    rho_steel = 7850;            % kg/m^3
    cost_conc  = Cc*(b*h);       % THB per meter length (unit depth 1 m)
    cost_steel = Cs*(As*rho_steel); % THB per meter (kg = volume * density)
    f = cost_conc + cost_steel;
end

Constraints

function g = nl_ineq_rc(x, fc, fy, Mu, phi)
    % g(x) <= 0
    b = x(1); d = x(2); rho = x(3);
    As = rho*b*d;
    a  = As*fy / (0.85*fc*b);
    Mn = As*fy*(d - a/2);       % N·m
    % capacity constraint: want phi*Mn >= Mu -> Mu - phi*Mn <= 0
    g1 = Mu - phi*Mn;
    % return column vector (could add more). Here one constraint
    g = [g1];
end

function x = repair_rc(x)
    % Example: enforce d >= 1.2*b as a practical proportion (optional)
    b = x(1); d = x(2); rho = x(3);
    if d < 1.2*b
        d = 1.2*b;
    end
    x = [b, d, rho];
end

Harmony Search

function out = hs_opt_general(prob, opts)
% HARMONY SEARCH (general, explicit constraints, feasibility rules)
% ---------------------------------------------------------------
% prob fields (explicit, easy to read/extend):
%   prob.obj(x)           : objective (minimize)
%   prob.lb, prob.ub      : 1xn bounds
%   prob.isInt            : 1xn logical mask for integer vars (optional)
%
%   % Linear constraints (optional, leave empty if none)
%   prob.A, prob.b        : A*x <= b
%   prob.Aeq, prob.beq    : Aeq*x = beq
%
%   % Nonlinear constraints (optional)
%   %   Return column vectors; sign convention: g(x) <= 0, h(x) = 0
%   prob.nl_ineq(x)       : returns [m1x1] vector g(x) (or [])
%   prob.nl_eq(x)         : returns [m2x1] vector h(x) (or [])
%
%   % Optional repair/projection function (discrete rules, spacing, etc.)
%   prob.repair(x)        : returns repaired x (or [] to skip)
%
% opts fields:
%   HMS, HMCR, PAR, bw, bwDamp, maxIter
%   seed (optional)
%
% Output:
%   out.best.x, out.best.f, out.best.feasible
%   out.hist.bestf, out.hist.bestvio
%   out.HM (final memory), out.rank (struct for debugging)
.
.
.
% ---------- main loop ----------
for it = 1:opts.maxIter
    xnew = zeros(1,n);
    for j = 1:n
        if rand < HMCR
            r = randi(HMS);
            xnew(j) = HM(r,j);
            if rand < PAR
                xnew(j) = xnew(j) + (2*rand-1)*bw(j);
            end
        else
            xnew(j) = prob.lb(j) + rand*(prob.ub(j)-prob.lb(j));
        end
        % immediate bound clamp helps stability
        xnew(j) = min(max(xnew(j), prob.lb(j)), prob.ub(j));
        if prob.isInt(j), xnew(j) = round(xnew(j)); end
    end

    s_new = eval_solution(xnew);

    % find worst (according to feasibility rules)
    iworst = 1;
    for k = 2:HMS
        if is_better(S(iworst), S(k))
            % keep iworst as is
        else
            iworst = k;
        end
    end

    % greedy replacement
    if is_better(s_new, S(iworst))
        HM(iworst,:) = s_new.x;
        S(iworst)    = s_new;

        % bubble-up to keep roughly sorted (cheap)
        k = iworst;
        while k>1 && is_better(S(k), S(k-1))
            % swap up
            t = S(k-1); S(k-1)=S(k); S(k)=t;
            tt = HM(k-1,:); HM(k-1,:)=HM(k,:); HM(k,:)=tt;
            k = k-1;
        end
    end

    % track best
    if is_better(S(1), best), best = S(1); end
    hist.bestf(it)   = best.f;
    hist.bestvio(it) = best.vio;

    % bandwidth damping
    bw = bw * opts.bwDamp;
end

out.best = best;
out.hist = hist;
out.HM   = HM;
out.rank = S;
end

Output & convergence plot

---

Reference:

A E Kayabekir et. al., Harmony Search Algorithm for Structural Engineering Problems, Springer, 2021